{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "# 第6章 逻辑斯谛回归与最大熵模型"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## 习题6.1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "&emsp;&emsp;确认逻辑斯谛分布属于指数分布族。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**解答：**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**解答思路：**  \n",
    "\n",
    "1. 列出逻辑斯谛分布的定义\n",
    "2. 列出指数分布族的定义\n",
    "3. 通过指数倾斜，证明逻辑斯谛分布的分布函数无法表示成指数分布族的分布函数形式"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**解答步骤：**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**第1步：逻辑斯谛分布的定义**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "&emsp;&emsp;根据书中第6章的逻辑斯谛分布的定义：\n",
    "\n",
    "> **定义6.1（逻辑斯谛分布）** 设$X$是连续随机变量，$X$服从逻辑斯谛分布是指$X$具有下列分布函数和密度函数：\n",
    "> \n",
    "> $$ F(x) = P(X \\leqslant x) = \\frac{1}{1 + \\text{e}^{-(x-\\mu)/\\gamma}} \\\\\n",
    "f(x) = F'(x) = \\frac{\\text{e}^{-(x-\\mu)/\\gamma}}{\\gamma(1+ \\text{e}^{-(x-\\mu)/\\gamma})^2} $$\n",
    "> \n",
    "> 式中，$\\mu$为位置参数，$\\gamma > 0$为形状参数。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**第2步：指数分布族的定义**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> 参考指数分布族的Wikipedia：https://en.wikipedia.org/wiki/Exponential_family  \n",
    "> 对于随机变量$x$，在给定参数$\\theta$下，其概率分别满足如下形式：\n",
    "> \n",
    "> $$ f(x|\\theta)=h(x)g(\\theta)\\exp(\\eta(\\theta)\\cdot T(x)) $$\n",
    "> \n",
    "> 称之为指数分布族。  \n",
    "> 其中：$g(\\theta)$表示归一化系数，$h(x)>0$ "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**注：**这里我们将通过替换$\\eta(\\theta)$来构造反例分布"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**第3步：证明逻辑斯谛分布的分布函数无法表示成指数分布族的分布函数形式**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "&emsp;&emsp;根据指数分布族的Wikipedia：https://en.wikipedia.org/wiki/Exponential_family\n",
    "> Other examples of distributions that are not exponential families are the F-distribution, Cauchy distribution, hypergeometric distribution and logistic distribution\n",
    "\n",
    "&emsp;&emsp;可知，逻辑斯谛分布不属于指数分布族"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "证明思路：   \n",
    "参考：https://stats.stackexchange.com/questions/275773/does-logistic-distribution-belongs-to-exponential-family \n",
    "\n",
    "1. 设$\\gamma=1$，可得单参数的逻辑斯谛分布；\n",
    "2. 计算当$\\mu=0$时，函数$f(x|\\mu=0)$；\n",
    "3. 根据逻辑斯谛分布的MGF（矩生成函数），可得$E(\\text{e}^{\\theta x})$；\n",
    "4. 根据指数倾斜的定义，证明单参数$\\theta$的指数倾斜密度函数无法表示成逻辑斯谛分布的密度函数形式；可证得，逻辑斯谛分布不属于指数分布族；"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "证明步骤："
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. 单参数逻辑斯谛分布：  \n",
    "设$\\gamma=1$，则单参数$\\mu$的逻辑斯谛分布：\n",
    "$$\n",
    "f(x|\\mu) = \\frac{\\text{e}^{-(x-\\mu)}}{(1+ \\text{e}^{-(x-\\mu)})^2}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2. 计算$f(x|\\mu=0)$  \n",
    "$$f(x|\\mu=0) = \\frac{\\text{e}^{-x}}{(1+ \\text{e}^{-x})^2}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "3. 逻辑斯谛分布的MGF矩生成函数"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "根据逻辑斯谛分布的Wikipedia：https://en.wikipedia.org/wiki/Logistic_distribution\n",
    ">  Logistic的MGF矩生成函数$M_X(\\theta)$：\n",
    "> $$ M_X(\\theta) = \\text{e}^{\\mu t}B(1-st, 1+st) $$ \n",
    "> \n",
    "> 其中$t \\in (-1/s, 1/s)$，$B$表示Beta函数。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "可知，当$\\mu=0, s=1$时，\n",
    "$$\n",
    "M_X(\\theta) = E(\\text{e}^{\\theta x}) = B(1-\\theta, 1+\\theta), \\quad \\theta \\in (-1, 1)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "4. 证明单参数$\\theta$的指数倾斜密度函数无法表示成逻辑斯谛分布的形式\n",
    "\n",
    "根据指数倾斜的Wikipedia：https://en.wikipedia.org/wiki/Exponential_tilting \n",
    "> &emsp;&emsp;给定一个随机变量$X$，其概率分布为$P$，概率密度为$f$和矩生成函数(MGF)为$M_X(\\theta) = E(e^{\\theta x})$，指数倾斜$P_{\\theta}$定义如下：\n",
    "> \n",
    "> $$ P_{\\theta}(X \\in dx) = \\frac{E[\\text{e}^{\\theta X} I(X \\in dx)]}{ M_X(\\theta)} = \\text{e}^{\\theta x - k(\\theta)}P(X \\in dx) $$\n",
    "> \n",
    "> &emsp;&emsp;其中，$k(\\theta)$表示为累积生成函数（CGF），即$\\log E(\\text{e}^{\\theta X})$，称$P_{\\theta}(X \\in dx)=f_{\\theta}(x)$为随机变量$X$的$\\theta$-tilted密度分布。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "&emsp;&emsp;综上，我们将使用$M_X(\\theta)$替换$\\eta(\\theta)$，可知\n",
    "$$\n",
    "f_{\\theta}(x)=\\text{e}^{\\theta x - k(\\theta)} f_0(x)\n",
    "$$\n",
    "\n",
    "&emsp;&emsp;其中，$k(\\theta)= \\log M_X(\\theta) = B(1-\\theta, 1+\\theta), \\quad \\theta \\in (-1, 1)$，根据指数倾斜性质，$f_{\\theta}(x)$无法表示逻辑斯谛分布的密度函数形式。  \n",
    "&emsp;&emsp;所以，逻辑斯谛分布不属于指数分布族。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 习题6.2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "&emsp;&emsp;写出逻辑斯谛回归模型学习的梯度下降算法。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**解答：**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**解答思路：**\n",
    "1. 写出逻辑斯谛回归模型\n",
    "2. 根据书中附录A梯度下降法，写出逻辑斯谛回归模型学习的梯度下降法\n",
    "3. 自编程实现逻辑斯谛回归模型学习的梯度下降法"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**解答步骤：**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**第1步：逻辑斯谛回归模型**  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "&emsp;&emsp;根据书中第6章逻辑斯谛回归模型的定义：\n",
    "\n",
    "> **定义6.2（逻辑斯谛回归模型）** 二项逻辑斯谛回归模型是如下的条件概率分别：\n",
    "> $$ P(Y=1 | x)=\\frac{\\exp (w \\cdot x+b)}{1+\\exp (w \\cdot x+b)} \\\\ \n",
    "P(Y=0 | x)=\\frac{1}{1+\\exp (w \\cdot x+b)} $$\n",
    "> \n",
    "> 这里，$x \\in R^n$是输入，$Y \\in \\{0, 1 \\}$是输出，$w \\in R^n$和$b \\in R^n$是参数，$w$称为权值向量，$b$称为偏置，$w \\cdot x$为$w$和$x$的内积。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**第2步：逻辑斯谛回归模型学习的梯度下降法**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "&emsp;&emsp;根据书中第6.1.3节的模型参数估计：\n",
    "> 逻辑斯谛回归模型的对数似然函数为\n",
    "> $$ L(w)=\\sum_{i=1}^N \\left[y_i (w \\cdot x_i)-\\log (1+\\exp (w \\cdot x_i))\\right] $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "&emsp;&emsp;将对数似然函数求偏导，可得\n",
    "$$\n",
    "\\frac{\\partial L(w)}{\\partial w^{(k)}}=\\sum_{i=1}^N\\left[x_i \\cdot y_i-\\frac{\\exp (w^{(k)} \\cdot x_i) \\cdot x_i}{1+\\exp (w^{(k)} \\cdot x_i)}\\right]\n",
    "$$  \n",
    "&emsp;&emsp;梯度函数为\n",
    "$$\n",
    "\\nabla L(w)=\\left[\\frac{\\partial L(w)}{\\partial w^{(0)}}, \\cdots, \\frac{\\partial L(w)}{\\partial w^{(n)}}\\right]\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "&emsp;&emsp;根据书中附录A 梯度下降法的算法：\n",
    "\n",
    "> **算法A.1（梯度下降法）**  \n",
    "> 输入：目标函数$f(x)$，梯度函数$g(x)= \\nabla f(x)$，计算精度$\\varepsilon$；  \n",
    "输出：$f(x)$的极小值$x^*$。  \n",
    "（1）取初始值$x^{(0)} \\in R^n$，置$k=0$。  \n",
    "（2）计算$f(x^{(k)})$。  \n",
    "（3）计算梯度$g_k=g(x^{(k)})$，当$\\|g_k\\| < \\varepsilon$时，停止迭代，令$x^* = x^{(k)}$；否则，令$p_k=-g(x^{(k)})$，求$\\lambda_k$，使\n",
    "> \n",
    "> $$ \\displaystyle f(x^{(k)}+\\lambda_k p_k) = \\min \\limits_{\\lambda \\geqslant 0}f(x^{(k)}+\\lambda p_k) $$\n",
    "> （4）置$x^{(k+1)}=x^{(k)}+\\lambda_k p_k$，计算$f(x^{(k+1)})$。当$\\|f(x^{(k+1)}) - f(x^{(k)})\\| < \\varepsilon$或 $\\|x^{(k+1)} - x^{(k)}\\| < \\varepsilon$时，停止迭代，令$x^* = x^{(k+1)}$。  \n",
    "（5）否则，置$k=k+1$，转步骤（3）。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "逻辑斯谛回归模型学习的梯度下降算法：  \n",
    "输入：目标函数$f(w)$，梯度函数$g(w) = \\nabla f(w)$，计算精度$\\varepsilon$  \n",
    "输出：$f(w)$的极大值$w^*$  \n",
    "（1）取初始值$w^{(0)} \\in R^n$，置$k=0$。  \n",
    "（2）计算$\\displaystyle f(w^{(k)})=\\sum_{i=1}^N \\left[y_i (w^{(k)} \\cdot x_i)-\\log (1+\\exp (w^{(k)} \\cdot x_i))\\right]$  \n",
    "（3）计算梯度$\\displaystyle g_k=g(w^{(k)})=\\sum_{i=1}^N\\left[x_i \\cdot y_i-\\frac{\\exp (w^{(k)} \\cdot x_i) \\cdot x_i}{1+\\exp (w^{(k)} \\cdot x_i)}\\right]$，当$\\|g_k\\| < \\varepsilon$时，停止迭代，令$w^* = w^{(k)}$；否则，令$p_k=-g(w^{(k)})$，求$\\lambda_k$，使\n",
    "\n",
    "$$\n",
    "\\displaystyle f(w^{(k)}+\\lambda_k p_k) = \\max_{\\lambda \\geqslant 0}f(w^{(k)}+\\lambda p_k)\n",
    "$$  \n",
    "（4）置$w^{(k+1)}=w^{(k)}+\\lambda_k p_k$，计算$f(w^{(k+1)})$。当$\\|f(w^{(k+1)}) - f(w^{(k)})\\| < \\varepsilon$或 $\\|w^{(k+1)} - w^{(k)}\\| < \\varepsilon$时，停止迭代，令$w^* = w^{(k+1)}$。  \n",
    "（5）否则，置$k=k+1$，转步骤（3）。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**第3步：自编程实现逻辑斯谛回归模型学习的梯度下降法**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.optimize import fminbound\n",
    "from pylab import mpl\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "# 图像显示中文\n",
    "mpl.rcParams['font.sans-serif'] = ['Microsoft YaHei']\n",
    "\n",
    "\n",
    "class MyLogisticRegression:\n",
    "    def __init__(self, max_iter=10000, distance=3, epsilon=1e-6):\n",
    "        \"\"\"\n",
    "        逻辑斯谛回归\n",
    "        :param max_iter: 最大迭代次数 \n",
    "        :param distance: 一维搜索的长度范围\n",
    "        :param epsilon: 迭代停止阈值\n",
    "        \"\"\"\n",
    "        self.max_iter = max_iter\n",
    "        self.epsilon = epsilon\n",
    "        # 权重\n",
    "        self.w = None\n",
    "        self.distance = distance\n",
    "        self._X = None\n",
    "        self._y = None\n",
    "\n",
    "    @staticmethod\n",
    "    def preprocessing(X):\n",
    "        \"\"\"将原始X末尾加上一列，该列数值全部为1\"\"\"\n",
    "        row = X.shape[0]\n",
    "        y = np.ones(row).reshape(row, 1)\n",
    "        return np.hstack((X, y))\n",
    "\n",
    "    @staticmethod\n",
    "    def sigmod(x):\n",
    "        return 1 / (1 + np.exp(-x))\n",
    "\n",
    "    def grad(self, w):\n",
    "        z = np.dot(self._X, w.T)\n",
    "        grad = self._X * (self._y - self.sigmod(z))\n",
    "        grad = grad.sum(axis=0)\n",
    "        return grad\n",
    "\n",
    "    def like_func(self, w):\n",
    "        z = np.dot(self._X, w.T)\n",
    "        f = self._y * z - np.log(1 + np.exp(z))\n",
    "        return np.sum(f)\n",
    "\n",
    "    def fit(self, data_x, data_y):\n",
    "        self._X = self.preprocessing(data_x)\n",
    "        self._y = data_y.T\n",
    "        # （1）取初始化w\n",
    "        w = np.array([[0] * self._X.shape[1]], dtype=np.float)\n",
    "        k = 0\n",
    "        # （2）计算f(w)\n",
    "        fw = self.like_func(w)\n",
    "        for _ in range(self.max_iter):\n",
    "            # 计算梯度g(w)\n",
    "            grad = self.grad(w)\n",
    "            # （3）当梯度g(w)的模长小于精度时，停止迭代\n",
    "            if (np.linalg.norm(grad, axis=0, keepdims=True) < self.epsilon).all():\n",
    "                self.w = w\n",
    "                break\n",
    "\n",
    "            # 梯度方向的一维函数\n",
    "            def f(x):\n",
    "                z = w - np.dot(x, grad)\n",
    "                return -self.like_func(z)\n",
    "\n",
    "            # （3）进行一维搜索，找到使得函数最大的lambda\n",
    "            _lambda = fminbound(f, -self.distance, self.distance)\n",
    "\n",
    "            # （4）设置w(k+1)\n",
    "            w1 = w - np.dot(_lambda, grad)\n",
    "            fw1 = self.like_func(w1)\n",
    "\n",
    "            # （4）当f(w(k+1))-f(w(k))的二范数小于精度，或w(k+1)-w(k)的二范数小于精度\n",
    "            if np.linalg.norm(fw1 - fw) < self.epsilon or \\\n",
    "                    (np.linalg.norm((w1 - w), axis=0, keepdims=True) < self.epsilon).all():\n",
    "                self.w = w1\n",
    "                break\n",
    "\n",
    "            # （5） 设置k=k+1\n",
    "            k += 1\n",
    "            w, fw = w1, fw1\n",
    "\n",
    "        self.grad_ = grad\n",
    "        self.n_iter_ = k\n",
    "        self.coef_ = self.w[0][:-1]\n",
    "        self.intercept_ = self.w[0][-1]\n",
    "\n",
    "    def predict(self, x):\n",
    "        p = self.sigmod(np.dot(self.preprocessing(x), self.w.T))\n",
    "        p[np.where(p > 0.5)] = 1\n",
    "        p[np.where(p < 0.5)] = 0\n",
    "        return p\n",
    "\n",
    "    def score(self, X, y):\n",
    "        y_c = self.predict(X)\n",
    "        # 计算准确率\n",
    "        error_rate = np.sum(np.abs(y_c - y.T)) / y_c.shape[0]\n",
    "        return 1 - error_rate\n",
    "\n",
    "    def draw(self, X, y):\n",
    "        # 分隔正负实例点\n",
    "        y = y[0]\n",
    "        X_po = X[np.where(y == 1)]\n",
    "        X_ne = X[np.where(y == 0)]\n",
    "        # 绘制数据集散点图\n",
    "        ax = plt.axes(projection='3d')\n",
    "        x_1 = X_po[0, :]\n",
    "        y_1 = X_po[1, :]\n",
    "        z_1 = X_po[2, :]\n",
    "        x_2 = X_ne[0, :]\n",
    "        y_2 = X_ne[1, :]\n",
    "        z_2 = X_ne[2, :]\n",
    "        ax.scatter(x_1, y_1, z_1, c=\"r\", label=\"正实例\")\n",
    "        ax.scatter(x_2, y_2, z_2, c=\"b\", label=\"负实例\")\n",
    "        ax.legend(loc='best')\n",
    "        # 绘制透明度为0.5的分隔超平面\n",
    "        x = np.linspace(-3, 3, 3)\n",
    "        y = np.linspace(-3, 3, 3)\n",
    "        x_3, y_3 = np.meshgrid(x, y)\n",
    "        a, b, c, d = self.w[0]\n",
    "        z_3 = -(a * x_3 + b * y_3 + d) / c\n",
    "        # 调节透明度\n",
    "        ax.plot_surface(x_3, y_3, z_3, alpha=0.5)\n",
    "        plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPMAAAD0CAYAAABO8xCHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAABu2UlEQVR4nO29d5Qc5bUtvqtzmpyTZiSNpNFIk6SRwBIgm2CyEBljY2PwfeB7DZdlY+DaLD/ez+H68vC1fW0D9rNsMMYGJBFMMGBMTpJGTNYETY49nXOsrvr90fpK1T0dqnuqZxR6r8UCZrqqq3tq13e+c/bZh2JZFllkkcWpD8lKX0AWWWQhDrJkziKL0wRZMmeRxWmCLJmzyOI0QZbMWWRxmiBL5iyyOE0gS/L7bN0qiywyD0qMk2RX5iyyOE2QJXMWWZwmyJI5iyxOE2TJnEUWpwmSJcCyyEIUBINBzMzMwOfzrfSlrBhUKhWqq6shl8szcn4qSaNFNpudhSgYHx9HTk4OioqKQFGiJG9PKbAsC7PZDKfTidWrV0f/OpvNzuLUgc/nO2OJDAAURaGoqCijkUmWzFksG85UIhNk+vNnyZxFFnGwsLAAi8US9/cMwyzj1SRHNgGWxRmB5557Dg8//HDC19x111342te+xv3/u+++i7feegv/7//9v0Wv7ejowJ/+9Cd89tln6O/vx6ZNm9DW1oZf/epXeOaZZ+Dz+XDrrbeK/TESIkvmLM4I3HDDDbjhhhsiftbQ0IDBwUHu/3/9619j06ZNyMnJ4X4WDAZx9tlnAwDcbje+/OUv44EHHsCBAwfwla98Bffccw9+/OMfY+/evXC73fi///f/oqamZnk+VBSyZM7i5IXRCExMAHV1QEnJkk9ns9lgtVpjZZMBAEqlEo888gjm5+exbt06nHvuuXjzzTdRUVGBpqYmvPfee+js7ARN03jmmWfw/e9/H/v378exY8dw77334qqrrkJ/f/+KkTm7Z87i5MRf/wrU1gIXXRT+91//uuRTfvrpp/jlL38Z9/f/8i//ggsvvBCXX345Hn30UfzqV7/Chx9+iPHxcQDArl27cM899+DJJ5+E2WwGALz44ov4xS9+gUAgkHB/vRzIkjmLkw9GI3D77YDXC9jt4X/ffnv45yJifn4e7e3t3D/z8/PYt28fbr31Vtxzzz345je/CafTGRF2A+G99O7duzEyMgKapqHT6TA0NIS2tjZRry9VZMPsLE4+TEwACkWYxARyefjnIoTbBBUVFejo6Ij42c0334zW1lZ8+umnOOusszA+Po5Vq1ZFvObRRx/FXXfdhTfffBMPPPAAvv3tb0Oj0SA3N1e0a0sHWTKvAFiWRTAYhEQigVQqPePrr4tQVwcEApE/CwbDP88gHnroIbz++uvc/z/yyCOYm5vDl7/8Ze5nW7ZswaOPPgoA+Nd//VfodDrU1dXhqquuyui1CUGWzMsMhmEQCAQilEBSqRRyuRwymSxLbiC8+u7dGw6t5fIwkffuFXVVjoWHHnoIDz30EADg6NGjuP3225GTk4MrrrgCd955J4qLixcds3fvXvh8Plx00UUAgCeeeAJ//vOfIZEs/w42S+ZlAsuyoGkaNE2DoihIpVLu5wzDwOv1ciT2+/3IycmBQqE4c8n9pS8BF14oWjZ79+7dmJubAwC0t7cDAHQ6HfffAFBTU4Pdu3fj+eefh9/vx6OPPorm5mY899xzuOyyy9Da2orvfOc72LBhA4Cw3vz111/H008/DQCwWCy46qqr4Pf7sXfv3iVdbzrINlosA1iWRSAQAMMwHDGDwWDc1/b19aGurg4ajQbA6bFyDwwMYOPGjSt9GQnx9NNPY3p6Gnv27EFDQ0PE71iWxUsvvQSVSoVLLrkk7feI8z2I8gfNrswZBk3TmJmZQSgUQlVVFSiKAsuyYFk2JikpiuJWbqlUGnPllslk3D+nKrlPRvD3xtGgKAp79uxZvotJA1kyZwj8sJphGC68ThWE3GQPxrIsQqEQaJrmXiOTybiVWyKRZMl9hiJL5gyAYRgEg0EurCarsVAkej05H0E0uSmKili5s+Q+c5Als4ggxCL7YbKaxiNnvFA7FcQiN03T6OjoQEtLS5bcZxCyZBYJpHYcCoUWESyazMlW61RX8uhjKYpCKBTi9tw0TXMPmCy5T19kySwCSO2YrLTR5FgKOZeKWCt3MBhcRG65XA6pVJol9ymMLJmXgOjacTyhQCwyu91uKBSKmOZumSQ/v8YNxCY3SabJZLKYD6dTHUK3NwsLC5DL5SgsLIz5e4ZhVkQcEg8nz5WcYiC1Y0LkRDcHn5yhUAhHjx7F0NAQuru7ceTIEYyNjcFqtSIUCkWcfznAL4ORMlcgEIDb7cbAwADMZjM8Hg+X0FupCEMsvP/++7j00ksFvfbdd9/F/fffH/N3HR0duOeee3DOOeegoKAA55xzDu666y4AwDPPPIMnnnhCrEsWjOzKnAZIkiteWB0NQma3242enh5UVVVh7dq1AMJ1aJvNBqPRiJGREchkMtA0DZfLBZ1Ot+xPfv7KbbVaUVpaikAgAL/fz0Uf0WH5qYKOjg58/etfR0VFBc455xwAQFdXF/7617/iyiuvxK9//Ws89thjWXOCMwFCw+poUBQFp9OJ7u5ubN68Gbm5udweWy6Xo6SkBCXH5Yp+vx9Hjx6FyWTCzMwMlEolCgoKUFBQAJ1Ot+whL2kGAU5EC4FAAIHjjRCZJLeY3gTPPvssvv/97+Oll17C5s2bAQC/+c1vsG3bNlx55ZUATn1zgiyZBSJW7VgIQqEQxsfH4fF4sGPHDshkib9ypVIJjUaD8vJy5OXlwev1wmq1YmpqCi6XCxqNhiO3RqPJKLmjQ2ryXstB7r/+NdxnoVCEG6j27g3LtdOFTqfDRRddhJtuugn5+fkRv9u5cyfkcjneffddBINBWCwW3HPPPejq6oLRaER7ezuampqwa9cu7Nq1C3v37l1kTvDEE0+suDlBlsxJEF07ToXITqcTfX19KCwshFqtTkpkAv4eW61WQ61Wo7KyEizLwuPxwGq1YmxsDB6PBzqdjiO3Wq1O70MmuZZkv+OTm+QSCLn5JTKhiSe+NwFpab799nDfRbor9OWXX47Dhw/jkUcewSWXXIKDBw+iv78ft912G3w+H6e33rdvH5566ik89NBD2Lp1K7773e8KNie499578cILL6R3gSIgS+YESFQ7Tnbc7OwspqensXnzZoRCIczOzi75eiiKglarhVarRXV1NViWhcvlgtVqxfDwMNdtVVBQIIoNbKrJrlhlMKPRCK/XC6lUuuh7jPedLoc3gdVqxczMzKKfZ80JTkMwDAOn04mhoSE0NzcLJjJN0+jv74dUKsX27dshlUpht9tFk3NGvy4nJwc5OTlYtWoVd81WqxU+nw+HDx9Gbm4ut3KnM+NoKWF89MrNbzLhvyaa3MvhTUCaX/jImhOcZuCH1aRMI/SGdjgcXPtiZWUl9/PlEo1IJBLk5eUhLy8PRqMRW7ZsgcPh4FYhhmGQl5eHgoIC5OfnJw37xbhmElqTc0V/l7HIXVREYe9e6rg3ASWaN8GePXtQUlKC2dlZ/OY3v8GaNWtwxRVX4Pvf/z5+9KMf4ZxzzsmaE5wuiO47lkgkgm5olmUxPT2N2dlZtLS0QKvVRvw+FpkzJefkQyqVcqsyEI4a7HY7rFYrJiYmQFEU8vPzUVBQgLy8vAgxCf9aloJoMic7PyH39dcDn/88MDlJYfVqCqWlFFg2/evR6/U4duwYfvazn2FwcBC//OUvcd5556Gnpwff/e53kZeXh4qKCrz//vtZc4JTHbEkmaFQCIcPH+bqjLEQDAbR19cHpVKJDRs2xCSEy+XC6OgoWlpaIo7jmxXwMTIygsLCwrjKIyE4fPgwtm3blvA1wWCQ85K22+2QyWTIz89HYWEhcnJy0NnZidbW1pifSShIriAUCi1q+BcCfsZcLpdzclN+aC4EL7zwAjo7O3HJJZdgx44di35/4MAB7u+YNSc4RZGodpxshbTZbOjv78fatWtRXl4e93VitkCKiegadyAQgNVqxfz8PIaGhuDz+TA9PY2ioqK0a9wsyy4p5CTvGQqFoFAo4u65k5H76quvxtVXXx33fa699lpce+21Ca8ja05wEiNZ7ThR6+LExAQWFhbQ1tbG2fvEQ6zznIyySIVCgbKyMpSVlQEADh06BIVCgenpabhcLqjV6pRr3GK0efIhZM+dzsp9OuCMJLPQ2nGsnwUCAfT29kKr1WL79u2CVp1oMpvNZvT390MikXB71oKCAi4htZJdVnxIJBJUVFRwNW4iYCEiGK1Wm7TGLRaZU91z839/ppD7jCNzdFidyh/XarXi6NGjWLduHUpLSwUfxy/JjI2NwWw2Y8uWLZBIJFy2eXJykktI+Xw+6HS6dD6eqIgmhUajgUajQVVVFac1t1qtOHbsGHw+H1fjLigogFKpBBDZWbRUYgs5Nhm5WZaN8E1bTnJn+gF9RpE5XUkmIaHJZMKWLVtSVlqRZNqRI0eQm5uL9vZ2hEIhMAwTkegiCSmLxYKRkRHMzs6isLAQBQUFyMnJWZFVJd57UhQFnU4HnU6HmpoaMAwDl8sFi8WCo0ePIhgMIi8vD8FgEEqlEgqFAhaLBYWFhcv6OaLfy+PxRGyLlmvlZlkWZrMZKpUqI+cHzhAyR9eOU0nIMAzDkXDbtm1pJXPsdjvsdjtaW1u5ZFO0YAE4kZByu92coojUiJ1OJ6fLJvLQkylklEgkyM3NRW5uLurq6sAwDOx2O8bGxjA1NQWpVAq3242FhYW0vsNAIACFQrHk6xRynkwRW6VSobq6WtRz8nHakzldSSYArpe3oaEhpmBAyHtPTExAr9dDp9NxRBYKlUqFiooKVFRUROiyR0ZGuFCcrNwkrD1ZIJFIOHEK+YeUwWw2m6AaNx9CSm1CcOjQIWzfvj3m74itMQFFnVouLKc1mZPZ+SQ6bnR0FDabDRqNJi0iB4NB9Pb2Qq1WY+vWrejs7BR8bDyRCV+XHR3W0jTNqbuISORkAPnupVIpioqKUFRUBODElsJkMmF0dDRC4JKbm7siCip+Lze59nguLCcjuU9LMqfbdwwAPp8PPT09KCwsRHt7Oz755JOU399ut6Ovr4+rP4dCIdGTH9FhbSgU4tRdk5OTcLvdGB0dRWFhIXJzc5ck/FgK4iW94tW49Xo9hoeHoVAoOHJHdy0tF2KRmxg1zM3NoaysDBqN5qSxWDrtyMyyLPR6PbRaLRQKRUpfrtFoxPDwMBoaGrgVJNX3npqawtzcXET9OdZKK5FI4nY2pVOakkqlEcm0Q4cOITc3FwaDgXMwWYlkmtAMdnSN2+fzcfkCl8vFCVgKCgqg1WrTuv54ijuh4JPbYrGgrKwswoWFrNwr5cJyWpGZWMrOzs6ipqZG8D6SYRgcO3YMTqcT27ZtW5QgEXJD0jSNvr4+yOVyrluKYCVEIxRFLXIwsVgsy55MS7ccFZ0vOHToECQSCSYmJuB2uxfVuIW8h5gGfKRPO55Rw7e+9S08+OCDaUlY08VpQebosJo0wwuB1+tFT08PSkpKsHXr1kU3BVlBE4Wp8bql0kUmRCNKpTJuMs3r9SInJycjyTSxRCMSiQRVVVWLatz8ZCAhd7zyj5hkjr4nots9LRZLRswiEuGUJ3Os2nGiEJaPhYUFjIyMoLGxMW7SKBGxWJbFzMwMZmZm0NzcHFfocTIlSYDYJgdOp3NRMi0YDIKmacEOKbEgBpmjSRhd4ybXb7VaMTg4iEAgENHqSSKtZA/lVJHoc7nd7mUX/pyyZI6WZPL/2FKpNCGZGYbB0NAQvF5vzLCaj3gPBpqmcfToUVAUtSisJjDbXegYnESuVo1QjHPEu9GXW85JUVTMZJrBYEB3dzcAcN1UqSbTxCBzsnPwr7+2thYMw0T0cYdCIeTn5yfV0Kd6TYlALJ2WE6ckmZPVjiUSSUxRBhB+Yvb29qK8vBwNDQ1Jb7RYZHa5XOjp6cGqVatiigBmTTYcHpjA2JwR5G9uNepRuXoD1laVwOv1ore3FzRNc6FtXl7eSWNbS5JpSqUSW7du5cpIfDtgocm0TKzMyUA07/n5+Vi9ejX3cFpYWIDdbkdHRwcXkgupcaeDYDAoisglFZxyZBZSO463ms7Pz2N8fBybNm1CXl6eoPeLXiXn5uYwMTGBpqamiJIJy7IYnTOhY2ACsybbovO4fEG8+EEXSnNUKFXR2NrSBKVSya2Ax44dg1KpRGFh4YrcCIkQyw5YaDJNrJV5KQ868nAi5aP6+nrYbDaYzea0a9wn29YJOIXInErtOJrMoVAIAwMDoGka27ZtS8kLi5yLf47t27dz+8hQiMHA5Dw6BidhdrgTXr/FYsHsrBdVlRXQzFrQVl+F0tJSrmnD6/XCYrHAbDbD7/fD4XBwJDmZFF7xkmmjo6Pwer0RyjSxVmYxyEMy0DKZDMXFxZwYKBAIwGazYWFhIWaNO1bzRiKI3fYpFKcEmaPtfJJ9Ufw9s8vlQm9vL6qqqlBTU5Pyl0xRFDweDzeJgpzDH6TRMzqDz4am4PL6E56DzE9mGCac7aYofNo/jr6xWXy+bT3WVoZXPLVajaqqKq4xIy8vLyIpReSPQvy7lgvJkmkOhwMTExMoLi6OaPNMBWJloeOdR6FQRDxU/X4/rFYrZmdn4XQ6oVKpOHJrtVruoZAIK0Hok+OOSIBUR8EA4dWU1JsnJycXhcSpwO/3Y3BwEM3NzcjLy4Pb68eR4Sn0jMzAH6STHu/z+WAwGCCRSBbJQh1uH17+qBd15UX4fNs65OtOiExiJaWItnliYoLTPhOLn5Nlvx193V1dXSgqKoLD4cDU1BRYlk15vyoWMYQ+FJRKJcrLy1FeXg6WZTkBy+TkJFwuF1QqFYLBIDweT8waN03TK6K4O2nJvBRJJsuymJubQ05OTkRInAoYhsHg4CC8Xm/Y+5qS4c1DR3F0ch6hUPKyF8uycDgccDgcqKiogF6vj/vaCb0ZT71hxdYNq7CtoZY7no9obTORP87NzcHpdHIuIGTferKAkJesetGabJlMFvFQikXaTK/MiUBR1KJBBGazGWNjYxHbCn6Nm4halhsnJZnT7TsGwlMkRkdHodVq0dzcnNb7k7C6vLwcrFyNNzqGoLe5IbRaxDAMjEYjAKCqqkrQDRRiGBwamMDApB6bqvJRWZi4rMGXPxIXENIH7fP54Pf7sbCwgIKCghVNpkWvqomSafGsiZaaACMIhUJLPg9FUVCpVNDpdGhsbATLnhhEMDQ0BL1ejyeffBI0TcNgMKRkYpHgPV8DMMey7DcSve6kIrNQO594xxIBx5o1a+B2x09GJYLBYMDw8DB0xRU4OGpE18Ao57AhBIFAAAsLC1yomWp46PT48FbnCMrytbiuoBgFOcnfl+8CQjqqDh06BI/Hg9nZWTAMw9WJM1WKiYdkIXK8ZBoZv5OTkwOlUinKhA6xRCP8PTNFRQ4i2LhxI8xmMx577DF86UtfQnNzM37+85+n/V4URV0MoBXAXLLXnjRkXkrfMdFFy2QybN++nQtvUwHDMBgcGkL/2CwcjBK2mbFUPwLXklhaWrpkR4lZkx1PvXkQW9atwvbGWihS2CqQyY2rV6/G6tWrubGx/NCWNGVkerJkKvvdeMm0mZkZ2Gw2HD58OELZleqEDoZhREkcJkqAKRQKNDQ0oK2tDX/84x+XJP6hKEoL4P8AeBjA5mSvPynITMJSspKlcnPZ7Xb09/dj9erVqKioAJC4IykWHC4XXnzzfczaA5CrNKAoH/c7IWosso8KBAKoqqoSbeVjGBYdQ5MYml7Auc1rsb6mLK3zRJdiSGhLJktqtVqulCT2fnspySuSTCMPx9raWq7Nc2pqCgBSMjgQa++dLFwn3yn5DEvALwH8NwBBYeGKkpkfVvf392Pnzp0pHTs5OQm9Xr9oioRQMnt8AbxzuBfvHu5FXkEh1OrFSYtkZKZpGgsLC1Cr1aioqBBllYs+h9Pjw2uf9qN3bA6fb1uPotylJVeiQ1u32w2LxcINn8vLy0NhYaFoklKxFGDRbZ6pJtPEJHOiB4cYUk6Kor4GgGVZ9jmKom4VcsyKkTnV2jEfgUAAfX19UKlUMe1uk2mzrU43OgYn8f6RPrjdHpSVlcUNvxKR2ev1wmg0ori4WNCeeqkllmmDFU//4xBa62tw9qa6lELveOA3LZDhc3a7HRaLBR6PB0eOHIkoJa1ECSweCWMl0/g90NHJNDGz4onILFI2+y4A+RRFDQLIA6CmKErCsuxt8Q5YETKna+cDnLC7ra+v5xrZoxFPmz1vtuPw4AQGJ+Zht5ggVyhRWVmZVMQfqxfZbrfD5XKhoqIi5b1bKBSCXq8HTdNQqVRQq9VQqVSLXC1igWFYfDY8haGpBZzbshYNq+JP0kgHpH5dUFAAi8WC5uZmWK3WRZLTpZgEpAqhD8Ho+jDJ8JNkGiGzVqtdkqIuFAol3HuL0THFsmw7+e/jK/M5J1U2e6m14/HxcRiNxqRTJKLD7PF5Ew4NTGDGYIXH44HJZEJxcTG0Gg1yNCo4vX7BISXDMDAYDJBKpZxaKxUQEQmRaAYCAXi9XtjtdgAQ3Gjv9vnx+sGj6DseehfnZaZDRy6XL5KcEuGK2+3mvLIzKTlNtz7Mz/CzLIuenh6u243vmZZqMi0UCiX8rG63e0lzwtLFspFZSO043hOYTJHQ6XSC7G4lEgnoUAhHJ+ZxeHACJpsLLMvCarXC6/WisrISMpkMLACHxwe5TAqFXAZ3DFkmf2X2+/0wGAzIz89PWVFGUVSEiIRsBVQqFScXDIVC8Pl8cDqd8Pv9CAQCnGAh3s02Y7ThL/84jJb6apy9aTWU8sz+SaMFFLFMBQsLC0WVnIqRhaaosNNmVVUVJ8l0OBxcIjAVZVqyBJjb7UZNTc2SrpcPlmWfAPBEstdlnMxCa8fxHD0sFgsGBgawfv16QVa1geOa6VcOH0NhiRNA+MtfWFiAQqGIGVYH6RCCdAhalRLBUAgBnkyToihuiLnNZkNpaWnKKxApu3m9Xk5EEmtPL5VKudCVuGd4vV6YzWZuNSBk4t9MDMui89j08ax3PTbWLj30FhKp8GustbW1XKuhxWLhJKd+vx92u31JktNMyDljjbu1Wq2CkmlC9swrMZEko2ROZRSMVCqNyBKyLIvR0VFYLBZs3bo1ad3W4wug89g0uo5NwxcIwu0LohAnklRFRUVJkxJunx8UBeRqVXB6fJziy+VycU/1VG9Imqah1+shkUhQUlLCHS/k5lQoFFAoFMjLy+M0wiQkJ0oktVoNpVIZbgjxBfDGoXDonc8rry0XorPNgUAAHR0dmJubg8PhgFqt5n6fiu/Ycsg5ZTJZ3GQav82zoKAgqfb6tCNzqpJMQmYgvK/s7e1Ffn4+2tvbE/4hbS4POgYn0T8+D5qX9CJhtdvtTilJxbLhBgiFXAow4XPI5XKUlZWlvDqQ/XlJSQksFovg4+K5j5BVGQhHG16vF06nE2azGTKZjPv9rMmGjrk5yPPK8LlNq6FUpJagA8RZDRUKBeRyOTZu3BhTcsr3HUskOV2ukhIfiZJpVqsVwWAQJSUlMXMFK+EyAmSAzOmOgiFkNplMGBoaSmp3u2Bx4NDABEZmDGCiQkLSckjTdFpJKgCw2Z0wmUwoKylCIEindA6WZWGz2eB2u7n9udhWQFKplCsp8cN4EpLTNI2Pe4YxMDmPXa3rsbG2fEUb6mNJTkmrZDLJ6XJ3TSW79u7ublRUVMDj8XC5gtzcXC4k54tG0gVFURIAbwCoBcACuJtl2TcSHSMqmZciyZRIJBgfH4fP50N7e3vcfenEvBmHBycwtRB7peO3HKY6DoZ8Bn6izO/3gw4xyNGq4HQnD12Xmu3mX4dQUBS1KCSfmZmBz+fDuM2GiakZrCovxiVnN2FNjTBhS6b7cSUSCfLy8pCXl5dUcipGgwQg7njZvLw8lJSUcLkCkky77777cPjwYTz66KO4+uqrcd5556VcuiRvA+CrLMvOUxR1CYAfI0zuuBCNzCzLwu/3p1U79nq9XDja3t4eU7kzNL2AjsFJGKzOuO8vtOUwHmIlyoiwxen2QSmXQSqVwOMLxDyeNFnk5eUhNzc34ndir8yJQCIisncNhUKwur344ysfoKpAje0Nq1BRVorCwsKMTiVMBYkkp2azGW63G16vd8mSUzEdSwj4ybQ//vGPuOSSS7Bjxw787W9/S9nZhoAN3yzzx/+3FkB3smNEIzMhcKpfFhEjFBYWLtqXBukQ+sbmcGRoEna3N+45yGpIvJXTeYr7fD4YjUYUFBRE7HciSlNBGggCOrUSviANmj6xR3e73TCbzSgrK1tyvVXsFZEfkvtYFh+PmrHRy6BoYSFuOWmlrG8I+JLTgYEBFBYWwu/3L5KcptNwsVQkC9d9Ph+uv/563HzzzUt6H4qi7gNwPwAjgIuTvV7UMDuVBgeGYTA8PAy3241t27ZhenqaS4B5/Scy015/MOF5EtV+hdyQ/BW9rKxsUSIm1vEurz8860mjgsPjhdlsgd/vT9hksZwrcyJQFIUQC/TN2lBemIvzWuqhkrAR5SSyb13q9Yr1eVmWhVarRVlZWYTkNLrhglx3piWnye6rZAqxFN7nYQAPUxR1DYA3KIrayCb4UldEzkma/8vKyrBhwwZQVHgKhc3pQf/sIPrG5iIy07FA2uPsdntcEib70oWaCMT6/hiGgdXphtViglqpSnmIOEkUxvujLwfx9RYH9r3TiU1rKnFO01rU19cjEAhwSSmXy4W+vj5u75pqSC7mHpV/Hr7kFDjRcMGXnJIacSYkp4nOl4m/G8uyz1MU9T8AigCY4r1u2cms1+sxOjqKTZs2IT8/HwBgsDrxXs84xvSWRXvNWBBCwmQrYTAYhF6vT2oiEO88xMmjsLAQWp0OOrUS/kAQwTiWQvzzkPcm4Ouzlzu0ZQH0jc1hZMaIHZvXoGlNJcrLy1FUVITe3l7U1dXBYrFgcHAQwWCQC28LCgqWzdQuWVgb3XARLTklbqFiGBwIQTrbzRjnWAPAw7KsnqKozwHwsSwbl8iAyGROFnoMDQ3B7/dj+/bt3D5n1mjF8+93wWyxCvqyEyWZoq8lHpnJ/laIiUCs8xA1WHl5ORcRuLx+SCUUcjQquDw+xHuM8DutFAoFJwYhbYhkuPdy3XgEvkAQb382hL7xOZy/ZQMKdSpIJJKIjqpohRdJ/BQVFcU0ORBzZU4ldI4lOTWbzfD5fJzBwVIkp0L620VAPoDXKYqSAjAAuDHZAcuyMrvdbvT09KCyshIbN26M+ANXlRTg1ks/hxfeOYihycQZ6FQklfG6nVI1EeCfh2VZmEwmrn4dfYOFGBZOjw8qhRwUhUX7fafTCZ/Ph8rKSq6uLpFIImyJgsEgp82enZ2NK+HMFAxWJ579ZwfW15SgSBa51YlWeBGV1PT0NJxOJ7cCEiGFmCvzUgwOcnJyoNFoYDab0draukhymorLabJrCQQCoiTkWJb9DMD6VI7JOJnJBIjNmzfHXUlzNCpcur0R1QUajJn9izLXDMPAZDKBYRjB2epoMqdrIkBeR47XaDQoLi5OeLwvEOQ+l9cfQCBIw+fzQS6Xo7KyEhKJJO7TWy6Xc40AxcXFiySchNipzp5OBSyA/vF5OO1WKPPL0Lw2dr08WiVFogsipMjJyeF0B0txXxFDARbP4CDa5VSlUiWUnC5TL3NayFiYTSZAhEIhQXa3UqkUFQU6nLt9Cw4eHcfhwQkwDJu2QR6fMKmaCESDpmnMzc2lfLzT4wNYBjazgQtJU7kpY0k4PR4P7HY7AoFAxKqdCZO+AM3gnc5h9I/P4/Nt61FZHH+kT7TJQSgUgtFohMViwWeffbYk3zExVvh4JEzmchotOU32YFopXTaQoZXZ6XSit7cXNTU1qK6uFvSHIGGnXCbFOc31aKyrwIF/foq+KX1aBnmk28lms6VtIkAy5sFgEDU1NSkfTx4iJSUlkCt8UAo8Pt73JZVKuS4l4tTi9XphMBjAsmxEIk1MGGxO7HvnCDbWlWPn5jXQqpOfXyqVIj8/HzqdDs3NzYt8x6JD8kQQY2UWoiITIjklUVM8Up82KzOREU5NTaU8RYLfaMEwDBZmp7C1tgBnNa3Hx31jcMdRXSWC2WyGUqlMS1bJz5iThoFUYLfb4XQ6uYeIye2GLBBAaa4OHp8fdGhpSRKKoqBUKqFUKpGfnw+GYeD1erlQl6ZprlMpTQUS/3+g1agwMmPE6KwJZ29ajZa1ybc7/BU12ncsug+a1Inz8/NjkkSMlTnVB0Isyen8/DzsdntEtMGfSXXakNloNMJqtaY1RYKQ2ev1cjVoMnJ1/apyfNgzgp7RGUFG9IFAAC6XCzk5OYtGwghBdNlqZmZG8LEsy8JoNIJlWW5/zP0O4dBbJpUiR6OAzZl40FwqIHY4Wq0WLMtidnaWS/iFQqGIVTuVm1qjDGfc+fO03us6xoXe1SX5CT9DvA6w6D5om83GdSXJZDIUFRVxdWIxkJZbickE6QcfAA4HmNpayD73OeTm5sLj8WDDhg2LpmHa7Xa88847S75WiqIUAH4B4EIAFIAHWJY9kOw4UclcWlqatl2KVCqF1+vFZ599hsbGRk4QAAAqhRwXtm/EptWV+OeRQSxY4ntik4y3TqdLK9yMLlulQirSu6zT6ZCXlxdxI/MTcnQoBKeHhkalAMueSJjxX7sUkDpnbm7uol5oq9UKiUQCtVoNjUYDuVwe8/1UChnUSjk8/tgRkcnuwv53P0PDqjKc21IPrWpxqCx0rxs9eoeQhMx28vv90Ov1KCgoSFsqmzKZ3W5IX30VrEoFFBVBMjICMAxCra1c5BAdbUxMTOCVV17BBx98gLa2Ntx///246aab0rncQgBvsyz7rxRFrQdwiKKov7Esm1AOuWx15kRgGAZjY2Pwer0477zz4va2VhTl4csXbUfXsRl82DsS4QgSXTay2+0pETGWrVAqnymdJJvHF4CEopCjUcLtDSxq5RQL0Yk0mqbh9Xphs9m4WdDk92qlAgq5DFaHC74AjWQSnsGpBYzNhUPv1vrqyEgkzcRVNEkOHjwIn8/HheSklJTKdI5UM+qU1QrQNHC8AsNWVkI6MYHQpk1xtwGrV6/Gueeei/r6enz3u9+Fy+US/H58sCyrB7D/+H8PUxRFI+ydbU903Iqb4Pt8PvT09KCoqIgruSQCRVFoW1+DdTWleLdzCENTCxErIikbpaKFDoVCMBgMXOko1RvQbrdz3Vrx9qfxriccevshl0qhUUjh8oZXwkzKOWUy2aJEmt/ng8tmgSEQdgxNZZ8doEN4v3sE/ePz+ELbelSXhqMqMRJXROpbV1fHTcPk2/vI5XIukZZIupnytcjloBgm/HegKMDvB6tUIsSygrLZCoVCFFM/iqK+DqCHZdmERAZWeGUmRgQbN25EYWEhFhYWBB+rUytxxY5mVBeMYP9bn3APA/61CCEEKX2lY9LH3x+n2q1F0zQ3vVEulyMYCiHoDUGjUiC+fkx8SCQUigvy4PGpEGJY6I6bCrpcLvh8Pq42r9FokuZBzA439r/XiQ2rynBuc31GOq+kUmlEq6TP5+MEIMQtlJCbvzCkSma2tBSh9eshHRoCe5y8oYsuEmTmJ8awOACgKOoBhJVflwl5/YqszCzLYmRkBDabLaERQbJzjI2NIeiy4j9uvxY9Y3ocHBjnxq3G887mg0zvS8ekL9H+OBb4DxcSkmu1Ws6ChiSoGIYByzLQKsMqskwt0BQArVoJf5AO18OPg5gKyuVy2Gw25Ofnc/3mQhNpQ1MLGJ8zoXFVMcq0mR1Sp1KpUFlZyUk3SSmpr68vwr2EpunUVmaKAnPuuWDr68OrckEBUFAAZno6qWe2GEk7iqJ+A0ALYCfLsh4hxyw7mf1+P3p6ejh/r3Se3MFgED09PdDpdNi6dSskEgk+t3kNGmrL8faRQUzozVydORZIljcYDHLSylRA3EzSEaHwS1b87YDf7+cEIUB4JdHm5EGnVadVlksEjUqBEMNEZKhjIdrBhGGYiESaVCqNsALm/y0DdAif9E9AQTHQFZahprQgwTvFR6qOK/xB73z3EqPRyJGwsLCQGxebEBIJ2KqqiB8J8cwWYTTN2QA2sCx7YSrHLWuYTWxzN2zYkFbJCAiToa+vL+ZEi4IcDa79/BYMTenx8gefwe9ffLOSaRJqtRrl5an7YjkcDtjt9rREKA6HAzKZjNuXk8iBoijOPxsIP/CMRiOMZjNm5wPIz9VCoVRDKpMvSemlVsgBiorrlMJHLBJF68hjJdI0Gk3EdA67x48D73ViXXUpzmupR45m+doo+e4larUa7PH97vj4ONxuN+fblcoM6+WYM4XwCNd2iqJGeD/7Fsuyryc6aFlWZpYNT6MwmUxJbXMT7W1mZmYwPT2N1tbWhKHMhlXlyDl/C/55uB8mL8tlicmKKsR2N9ZnIKFmqvtj4hGlUChQWlqadD8vlUohk8lQWlrKW7W9AOOEN0CnrM9WyGRQyKVJV+JUEZ1I8/v9EdM5+EaGx2YMmNCbsX1jHbasr4FU4Pcn1r47FApBo9GgtLQ0IiQ3m82cuoufJY/39xWyZ15qmM2y7OMAHk/1ONHJHH2j8qdRJLPNJVMeFnUjhUI4evQoWJbF9u3bBa1OapUS7esqUVhWhX8eGcTQ+HTaKyrLspibmxPUZBEN0vfM97dOhuj6NH/VLqKAYMAPkzW5PlsulUCtVMDl9SNA00gVqXxO/nWS6Rx2ux1erxezs7Nc+eu9ziEu611bnjzbmynPbH5IzjcUNBqNGBkZ4bLR0SH5GafNJrDZbOjv78e6desEZfiICoyfYPB4POju7kZ1dbVgnTdwwsKoJF+HpnI1FHQh5gsLEKATJ8Wi4fP5EAwG01rNSYKtvLwcXm98D7NUEGIBiVyJNbU18AeDcLm98Hg83HB5tVoNrVYDtVKGEMvCKfJqLBRSqZR7eOXn53NWwCaTCQsLCxidnMamNdW45HPNyM+J/72m2sscD8keCtGGgnyfbK/Xy2XJkxngE+XhSiAjZGbZE7OTt2zZIthNka/PBk6Y/W3evBl5efE7dmJBIpEgEAjg0KFDqKiowLWXbIbXH8R7XcM4OjGf/AQ4sT+Wy+UpJbqIAIXfu+zz+VIyHEiW+HF5/ZBQ4bKS87hNDsMwkLA0LDYnHC4vaMbAlZVS3WuLVecmSb5YibT+sRl0D41hc20pdjSvR2lJ8aLE1FJ6mflIdYVXq9WoqqpCVVUV5xNnsVi4KaRk1Y4OyT0ez+lDZpqm0dPTE3d2ciIQMjMMg5GRETidTmzbtk1wcoIPh8MBs9mMrVu3ctJQjUqBS8/ejM1rKvHPjkGYHbG10dH747m5OcF7N4ZhsLCwALlcLtrw9bjvxYbNEJRyGVRyGbzBIAJBCQqLiuD1+ZCXlxfRVUXC8XTC/XQQ74EQnUibcwXx0icD2FiZiwK1jEtMEauf5Z5mEQ2KoriGC5vNhoaGBrhcrkVjbkOhENxud1pttlHvpwRwG4Avsix7tdDjRCfzwMAAKioqUF6e+vAysoINDg6ioKAAW7ZsSfmGIhpZvV7PjeyMRk1pIW65+Gx0DE3i0/7xCPNAku3m74+FClBIg0YsAUoqirRUPrNGKQcLwO7xQadWgmWA4PHPo1AoFnVVOZ3hSR18CadY0xpjQchnITmMAYMPdeVFKM/Jh8vl4qYzMgwDh8OxaIBbKhDzoaBQKBZ5jlksFnzve9/D9PQ07rzzTlx88cW47rrr0r3eIQCdAFJa4kX/KzY1NaUdogWDQRw9ehSbNm1Kq3RF0zT6+vqgUCjQ1taGvr6+uK+VSiU4q3E1V5semzPFzXYLISIRgqTTe50OVAo5pBIJ3L4Te+KwBXDYhyz6ForuqgoGg/B4PDAajWAYhiM2MRUUI8xO5xwTejOmDVZsa6hFe1sbHHY7JiYmuM4krVbLdVQJEvq4XJDMzIByODK29yYh+Z/+9Cece+65+PrXv46Ojo6lRDatx/95MJWDRCezEOVVNMge22azob6+Pi0iE5+x2tpaVFZWgqZpQdeRp1Xj6vPa8PFnfXjlw84Ikz6CZDc3EYLwGzSikSpB4toKyaRQKWRwemIntpjjPmRymRQalSKm7zh/DxurF1omk0Emk4liKpjODR1iGHx6dBwDk3psra+AVqvFhg0bOGsis9m8qOkiPz9/Ecmkb70F9S23ABSF9mAQ7l/9Ckiviyni8yT6TBKJBDt37sTOnTvTfg+WZW3pfG8r3mjBX02rq6vT2teQvUtTUxPnMybUkJ+Y8aupIB78l+vRMTSFI0NTgjqYEvUuJwORHZLSEj9JFesPKZVQ0KjCZSanJ/lDKkCH4PEFkKNRwRcIgo5jAQzEXrUdDgd8Ph9mZ2cXrdpCsdQasd3txWsHjyJfJUVZVQ3ydRrOmqi2tjZC4TUyMgKlUsmt2mqahvorXwHlCSshpQBy/v3f4d61C2xFRdrXlAgsyy6L33k8rCiZXS4Xenp6UFdXh8rKSkxOTqa0qhONt91uX5QoE7ISBgIBdHd3o7CwkDPjP691PTbWVeCfHYOYNdninovsrbVabUrabL6UtKqqiivZRCepyPtRAHQaFTz+QNzVOBGEWgDzr1OhUECr1UIikSA/Px8+nw8ej4dbtfkSzkQQQ/DBsizmLE489cYhtG9YhfaGWshl4YderHKS2WzGyMgIpIOD2IYwiblzyeWQjIwgtAQyJ7qnSN1/pZAR0YgQzM/PY3x8PMJeKLo0lQhEjJKbm4utW7cuet9k1+FwONDb2xuzBl6Sn4MbL2hH//gc3u8+tojMRAiSau2ZZVno9XooFAqUl5dzOl+SpOLPXA4EAnA7bJAplKCX6G4ZYQEMwBtIPPKHXCuwOPNM9tpCHUzEIDNFUQgxDA4OTGBgUo/zWupRX71Yt6BWqzk9AlNdDSo6MvP74SwpgToD3VzAyvp/ASuwMjMMg8HBQc4Mn7/HlEqlCAaT32gOh4PTZ6fTbkYeJIlkoRRFYfOaKqytKsGf//ZPGI+HtnwhSCols1AoBKfTieLiYk7+GA0y4K20uACTk9OQKdXweDzwWqwAAI1GsySbXc4CWK2Exx9EKI09sVwu58o08RoviINJJpJoDo8Pr3zSh9ryQny+dT0KcmKXgSSFhfD/5jdQfetbgFwOxueD+b77MB4KwXXoEHQ6XdjdRCaDXCoF8vPDfcsJkCwjLsZc5qVgWcns8/nQ3d2NsrKyRWb4QPhmppPIDufm5jA5OYmWlpa09NVkWJ1QnzK1UoGdm2rByjX4+yfdcDgcKXdakdKFRqNJKChQKeWQHG+ECNAnVm0ijyS6ZxLOEXKnmqV1ev2QSiXI0SjDM6fj3MTJHhixVm3yWWma5iSeSykLxQvVJ/UWPPXmQWxdvwrbN9ZxoTcf9PXXw71zJyQjI+iy29F05ZXYdPycLrsd9F/+Au+nn8LNssDmzWC/9jXklZYm1GUvl5STZdl3AbybyjHLFmabzWYMDg5yRgSxkChpRVb0QCCAbdu2pVwbDQQCXOtlW1tbyokcl2UBl2/bABdU+KR/DEGBslAyYbKoqAg+X+SwdvJZVUo5FDJZwlZH/ljW6KYGYguUyNMrGqEQA6fHD3UcH7J0IJfLIZfLkZubywlvgsEg5ufnE7ZLJkN89xAWhwcnMTi1gPNa6rEuRujNVlYiVFkJ3+HDEefLGxqC/NgxMFu2gGEYBEdGYHjlFXQ0N3NG+NGGFyezAT6wDCszMRGwWCxJjQji7ZnJil5aWhpzRU8Gp9OJnp4ewRpxPjweDxYWFlBVVYUNGzYAADasKsM7ncMYmTHEPY4kumiaRmVlZUQ7JhFCKBQyqORyOD0++I6XkJKVPshr+E0N0a2ISqUSDMMkvfmAEyN0wqF3ACHmxCiepYCiKMhkMi6CINdotVpB07TgsTtCkmhOjw+vftKHmtICfKFtPQpzkxOKmpkBq9EAEgkkEgmUFRWoCQZRvn07l+zjz4IuKiqCUqlM2jG1Uk0WQIbJHMtEIBFikZloYRsaGjj3xlRA9sctLS0pf9EkmigpKYmIJnK1alx1TgtGZ414+7NBONyLV1y9Xg+VShUxQJ6QWHL8HB5fAB5/kPteSLY7GAxyrxVC7litiB6PB3q9nguFk62ITq8fMqkEORpFOPQWAfwHQvQ1xhq7EyuySCUjPm2w4s//OIS2dTU4q7EOigTRG1teDsrj4Ty+KLsdoZYWAFhkhG+322E2m7kqxOTkZMzJHKfdnpl8uEQmAvHAJzPLspiamsL8/HxKzRoE5KaenZ3Ftm3bUmp7JO+t1+vR3t6O6enpmOH/2qoSrCorxKdHx9FxfJwOkXQWFBQsfniwLHI0Snj9wQirHvKdURSFQCAAo9HIPTz470uInejmJqu2VCpFZWUlN9KGrIiJss80L/Smg4FFyUiGAdxuQCIBNJqk+aKIzxb9s3huodH5gFTLWwzD4sjQFIamFnBucz02rCqLGWUw27cjNDICaU9P2CKorg6hCxcbe/BnQRcXF2N+fh4KhYKbzMH3HBNrZaYo6gYA/wUgBOAnLMv+QchxGVmZhZoIRIOQORQKob+/HxKJJOVmDeBEREBRFJqbm1MiMsMwXO/0tm3bIDkehiVSZJ3bXI/G2nK89N4RdA9OLZJ0siyL/BwtHDYrhkbGuSd/dDac1ElLS0u535H3JSs1V38WQGwgvCKSvl2yIhJyS6XSiFWbuw5/EP4gDa0qnJBjWBZ+P9DZKYHLRYFlgcpKBo2NbEJCCw3VE5kckPIdUawJJbbL68ffD/ajb3wO5zavXXwPyeWgv/pVhAwGgGXBFhcDSfIwRJfNtwEmBgePP/44nnzySaxZswbnnnsuzj777LQ07xRF5QD4GYCzESZzF0VRL7Msa0x2rOhk9vl8sNvtgk0E+JBKpVzbYnV1NWpqalJ+fzLnau3atZiZmUlp7+f3+9HV1YWysjLU1tZGrJjJ1GQehxWNpUq0bbgAnx6d5Mzj1Qo5GJaBxxdASVl5xEpJjPw0Gg2CwSDcbjcqKioivjf+NZCHCvknnVWbvyKS7DO/Zkxsf1iWhdtHQ6OVQKOQoq+PhttNITeXBcsCs7MSFBczKCtL/P2mmt+IzgdYLBawLJv2sLxpgxVP/+MQ8qQ0NjfTUMp5tzxFgRUYNQKJDQ7uvvtuzrz/qaeegtvtxsUXX5zSZz+OiwG8x7Ls7PH3eBvABQCeSXag6GRWqVTYtGlTWsdarVY4nU6cddZZKfcvA4Ber8fY2Biam5uh0+kwNzcnWIRCtgWx/MmSZdnJEPlt27ZBKpWioa4KH/eNYGTGCLsrLCckN3X08Defz8clypRKJddCl0jjzd+Dk2vgSwljEWhqioLBQKGwkMXq1eEVlZ99JjVjkvwBwium1+9HMCSDz62GVh0CQIOiAKkU8CTxjBSjzkzIrdFowLJhj+9oMwYS5cR7cIRCDIYXLPjT6wdxTvNabKxNvaMvfJ7EpSmGYbBr1y5cf/31aZ3/OGoATPL+fwaAIMnaiinA+OBnvDUaTcpEZlkWx44d4/qfScgoVJ+dTEQSTxpKwvm8vDxuLhYQHu3y+db1WF9dirc/G4LB6ox73Xa7HVqtFgUFBTE7mTQaTdz+Y/4DgpyP3PAkmw0AH3wgR0eHDCwb3u9u2hTCJZdEPuSia8Y2mw1+v5/r62Yl+bBZtSgv0SBABxAKsRCyPRRLAUbOFav27nA44Pf7OUPB6FWbuJW4fX68cego+sbm8Pm29SjJT21/m4zMZLrlEqEAwL9pGYTD7aTIyJ45lQ6hYDCI3t5eaLVabN26FZ9++mlK70UIlZOTs6j/ORmZyUPA5XIlFJHE6gTzeDzo6urCmjVrInq3SfjLMAwqivJw84Xb0D0yi4/6RiPG6ZDh7Xl5edwNwFdXsSwLj8fD9R8Tt5NEriEURXEELCsrg1QqhdPJ4tAhGfiX39srxbZtDIqK4v+NiBkfUXrpdD50dXkxMe2FTCpH/RoJCgoYJLqFxNJmx0N07T16xC0Jx6O3H7MmG/761mG01Ffj7MY6KBXCciqhUChh/kUkZ855AJ/n/X81gINCDlzRRguyv40mhFCQRo14xxODwFggjig6nS6piEQikUQo00jJKtrOiE9k/h72xDidYQxPL3BWusXFxXF7nymK4jqZAHDhJZn6wV+1+d+H3W5HeXk592Dy+yWQShFBZqmUhcvFID//xB4w+vPzSRQebarBeecBPh+LUIhGKOSEyx5AkKYhlStjRhBiyTmFuqJE69x9Ph+cTicn1iETRGQyGRiWxcDkPKqK82LqvGOB5BXiQSTRyBsA/pOiqFIAEgA7ANwh5MAVI3OsRotUsLCwgNHR0YTHx1uZ3W43uru7sXr1alQI6KDhJ8Cmp6cxNze3yDKYZdlwOHr85ou+AcPjdJpwpE+Bv71/JIJwQsDvP44l7STvX1FREZGkKShgIZMBgQhxGYXS0sWJvegk2uLmFUCtpgDIAeRxn1vKhmC02CIcTEi4LmaYnQrIZA6tVss1r4RCIa5ldWNtBS4+qwHlpcJ754UowJbq/8Wy7AJFUd8H8MnxH32HZdn4s395WPYwm/QPezyeuPXfRH/A6LbHRGFPrPCYzLdKxSSQnOfo0aMIBoNob29ftCcjK3e8MhqxM6ICbjxw2zXoGpnD4cGJhH3G8cAPLxmGgcFg4LTQCwsLXDgul8shkwFf/jKNffulsNso6HQsrr02BI3mhFCF/JufIRfam0tRFBhKhoryMsikUticLng8HhgMBgQCAe5a020OIe+xVMhkMuTn56OmshznNa9FvkoCs9mMmakJqNVqbqRsIoXiMu2ZwbLsEwCeSPW4ZV2ZyWiawsLCuKEteRDE+h0Jjcn+WkgjAP/mnJycxMLCQsrzrUKhEObm5lBdXb1ITkoMCBOVhRiGwcDAACQSCVpbWyPH6Xw2hEm9WfC1RJ/XYDBApVIhPz8fQPg7Ii2KNB02zNfpNPjmnSoA1KK6cPQqTB5MLpeLiwLI7xN9xgAdQoAOoTAvBzqNGgE63O8tl8uTJqgSQax9t4Si0FpfjR1NazhlWElJCZebiHYvKSoqWuS8KWSaxUo5cwLLSGbiob1+/XrOCC0WiHAkeoUj+2OhoTE5F9nD9vf3g6IoTggiFG63GyMjI8jJycHatWu5n8faH8cC6bsuKSlBTU1NxOsKcjS4dlcbhqb0eK/rWEoTJ2Il0IDFIhFiB2Q2mwUl0UKhENerTbYRsUpf8T4zsQDO0aigp8LSSHIt0WUlIS2dYpA5T6vC9tY12NayftHv+LmJVatWcSNjiXuNSqXi3EuSTbPwer0pKxXFxLKQeXp6GjMzM2hra0tqQ0rIzA+fyRfb3Nyc0pNPIpHA7/fj8OHDKC8vx6pVq1K6MUiia/Xq1RGDswmRyR833jndbjcnYEn0ANuwqhx1FcX4uHcU3aMzYJjkDilkcF2ihAxFURHlJkImg8HAlb60Wi0XpZDEXElJSUTkEl36SiZYIRbAcqkEWrUC/iATt6yUrKVzKUk0qUSC7RvrUJUnByuwd5s/MpY8DM1mM4aGhmC32yGVSlFWVhZz0Dspga0UMrZnBsJP+YGBATAMI1gRFq3PHh0dhc1mS8s/m2izm5ubU27SILrw9vZ2rkRErokkuhIR2WKxcHtzIQ8gpVyGL2zZgMbVFfhnxxD0FnvM18WSfApFLBM/EgITY4iysrK4WxChghXyTyAYgscXRJ5OA3+A5iyAAeEtnemuzFXF+bhg6wYU5moxOzsLKg2S8R+GNTU1+Oyzz1BYWBjbc+x4CUxMUBR1FoDfANjDsuxMstdnbGX2er3o7u5GRUVFSisiITNN0+jt7YVarcaWLVtSfuLNzc1hfn4eVVVVKRGZ9E3TNM0lusg0imQZa4LZ2VnMzc1hy5YtKXtClRXk4ksXtqNndBYf9o7Cz+szdjqdcDgciySf6YBv4ud0OjnxislkiriJ4z0w4glWYiXRnB7fcTOE2D5k0RJOInklwhUiJhJixKCUy3BO81psXl15IlJgGFG8wRmGQUlJCdc4RB6sw8PD+Pd//3e4XC68/vrr2LVrlxhG+A8jbLdbK/SYjJDZbDajv78fjY2NMU3oE0EqlcLtdqOvr48z+ksFfDeRtWvXIhAQPts4GAyiq6sLRUVFWL16NXczkEQayRgnyliPjIzA4/Fgy5YtS5qg0FJfjfqqErzffQwDk3pYrVb4/X5uHKxYsNls8Hq9Ee6i8fTjiVaf6FXb4XBwPyO5BTsdgkohh0wq5bTrscCXvM7OzkKn00Ws2vFaOuurS/D51vXQqSMfoMlmKgtFdJTA9xx77bXX8IUvfAFvvfUW9u3bhz/8QVCjUyL8hA1b7k4IPSAjZJbJZElHt8aD3+/H8PAw2traONtcoSBqsNzcXLS1tcFoNC5y94gHkmBbu3ZtRMsmCaetVitGRkZQUlKC/Pz8RTd1KBRCX18ftFotmpubRSGcVq3ExdsboaDd6PC6Un4wJgNpsIieUx1LP8535ySrdrzVjrQy8hOVZJX2HiexTq1EgA5xZgiJQNo2AcQ0OSgpzMeln2vG+lWxhUdiTbMA4pfJKCo8IO9nP/uZKO/Dsqwt1WMyQub8/PykXl7RIPpsh8OB+vr6lIlMhCB8NZhQQ35Se45OsJ1wBFHg7LPPhtVqxfz8PAYHB5GTk4OSkhIUFRVxJbPq6uqUI4lEIFLX+ppy7Dq7HUeGp3BoYEKwZVE8sGzY75vMgE6EWJ1WifTjZrMZDMNwc6j55yHvHe7ICkBCUdCo5Mez+PEz5Pyf8dslwbJYU56P+lId7AvT6LLouXoxP8wVk8zxkK7664477sDvfve7jqgf386ybHeq51pxE3wA3P5YpVItKt8IASEj3wQfENZoQaZVRteeo/fHMpmMmy9EwkiTycSN/KysrBR15SRWSXV1dVykcFbjamxYVY53PhvC+LwprfOSwXbpNLQAifXj5MFXUlIiOBz3+GmolUqwLMNZGCWraQNAYa4WF2zdgKrifO5nfN9sn8/H1YuTlZTEQLouI7/97W/x29/+tl2Ma8hoNlsIyPxlMlZmampKcNsiUVUZjcaYQpBkrYsDAwMIhUKLas/JEl2kAYG4grS2tsLtdnOGg4WFhSgpKRFkjh8LTqcTfX192LhxIycGIcjXqXH1ea0Ynl7Ae13HIhxLkoEY9+fn54tib0NqtBqNhpt8KZVKE+rHo48HwoITigIKcrRwef2gjwtxgBPREfk7SCQUtjfUYdvGWkhjzHwie9hQKMRNu1hYWIDL5UJ5eTmXeU4VyUpkK23mB6zwymw0GjE8PBwhrRTqnU3cSKRSKdrb22M+eeM1WpBJFsXFxairq4ur6EqU6JqamoLJZMLWrVshl8tRWFiImpoahEIhmM1mzM7OYmBgICIcF5JRJWWP5ubmhDfH+poy1JUX4ZP+MXQem05amw4Gg5wYRExhA1npdTodt0WJpx9PlJFm2RM+ZDqtCg6nByazmbMOYlkWFUV5uGDrBpTk5wjykyMhN03TKC0thdfrxdDQEILBYFyVV6LPmejBvNJmfsAKkZllWYyPj8NkMi1aUUkpKBFICErKXvEQa89MEl3RBvpCFV2kdMWyLNra2hYPK5NKUVpaitLSUi4cNxqNmJiYgFwuR0lJCYqLi2MSamZmhvM8E1JDVshl2NW6Ho11FfjnkSHMHR+nEw0iBkmnNp0IiVb6WHVkUm7i905Ha+vpEAOHywuH3QatSglNTi4Uchl2bl6DTXXl3PuSyEkqlSZ88JLXa7VaFBcXcw9cvspLiDZbSJNFJsjMsmyd0Ncue5hNBsUplcqYK2qitkXghCNIIv9tgugwm0QC8RJdyYhMElKFhYURtkLxQMLxvLw81NfXw+v1wmQyYWBgAMFgEEVFRSgpKUFOTg7GxsbgdrvTKmmV5OfgxvO3om98Dh90j0R4YJMSU6pdWslA0zT0er2glZ5fRybHRuvHiV0REFb8KZVKaHNzsbG2Ajub1iJHc6Iywq/5k78vTdOcX1v0PRWdAItWeUVrs4lndm5ubkS57mT2zAaWeWUm++NVq1ahqqoq5muie4f5INMshMhCybnIH35ychIGg2GRkkyoEMTj8aC3txerV69OayQOEN7T1dTUoKamBjRNw2w2Y3p6GkajEUqlMkL7nSooikLTmirUV5Xgg+4R9E/Mw+FwwOl0LmqLXCqInDRa9ikUifTjoVAIarUa5SVFuHBbI9bXLPboIp+FkItIa8m/STRGFHqJstnR2myapmG1WqHX6zE0NMTNg07WHCJWx9RSsGxkFtp6GGtlJkIQ0jYpdIUharL+/n6wLLsoEhBKZJvNhoGBAWzatCnlklk8yGQyFBYWYmZmBmvXrkVubu6icLykpCTlWr1aqcBF2zYiR0rjnU4rKioqRBWZ+Hw+zskkFdfTeCAiEJVKBb1ejxydDuurS7C6SAXb/CRG/E6UlJQgNzc37ufgr8YkwiL/sGykD3myh1p01YI8ZKanp+H3+zE2Nsat2vzrOSP2zMkyztGINsLn+2y1tramdGPSNA2Hw4HS0tJFiS4hPchA2ERhenoabW1taYlg4oHIXdesWcOt9CRz7fV6YTQa0d/fj1AohKKiIhQXFye8oQlYlsXg4CDy1DLcd+sedI7M4mD/OAIp1v1jgQhHxA7ZSafWqsoyXHP+dlSXhEt8wWAQFosF09PTcDqdyMnJQXFxMYqKiuI+SPjEDoVCGB8fh1arBUVRnEyYoihIpdKkxKYoitv35+XlQa/Xc0aRg4OD3PC5wsJCeDyepNu+VEBR1NUA/j8AOQD+CeB/sSybsMyT8UaLvr4+yOXyuBnnaPDJnMwWKBFcLhe6u7uhUCiwevVq7udC98ekycPlcmHLli2i3rx2ux1Hjx5FY2NjzChFrVZj1apVXNhHVgan04m8vDxuwkZ02Ee+75ycHE6Ouq2hFhtqyvBu5xBGZpNaL8cFsSQSQxcefc2GhQV8oX0TLt7RChnv3HK5HGVlZSgrK4uo7U9NTUEikXD7XkJWPkiS1e/3o6mpKUJWSkLx6HA8WRJNLpdHJDddLhdMJhP+67/+Cy+99BK2bNmCpqYmrmd9iSgAcBaAIIC/A7gRwF8SHUAlqZ+l3X9mt9vR2dmJmpoaVFdXCz7O7/ejt7cXtbW1MZNVQsBPdPX29mLHjh0AhBOZhOYqlQrr1q0TNUw1GAwYGxtDS0tLWlM6SO3UbDaHpYzHs+NSqRTd3d0oLy+P+32PzhrxTucwHG5vSu/rcDjgdrtRVlYm6t6bpmkE3Xbc8MUdWFeXmkc6MS40mUzweDzIz89HcXExCgsLIZFIMDIygmAwGHc2GQnBCbkJCKmjP6fBYIDH40FdXV3M67n33ntRUFCA2dlZ3HvvvWg5PupGIBLeYBRF/QzAFMuyv0z0uoyszCzLore3N6boIRkkEgncbjcmJiZSbntMlugSQmTihlJRUZHSQ0gIpqenYTAYuNp0qqAoihuVsm7dOk5W2dvbC6fzxN4yXtsgf5zOZ0NTguYzW61WBAKBRfrtJYNlUKFhceXlF6d8jwCAUqlEVVUVqqqqwDAMbDYbjEYjRkdHueaQzZs3J9xnA+FIUC6Xc6QmeZTo0leybHYoFMKVV16JnTt3pvxZEoGiqEoA1wD4YrLXZizM3rZtW8rHEZ8tmqYFDZrjg4yVAZB2osvlcqGvrw/r1q1La0hdPLBs2NLX7/fHrE2nC41Gg6KiIszNzaG5uRk0TWNychIulytuOM4fp/PPI0OYMVrjnp/orIXOChOKysIclMj92L61TRSbHYlEgsLCQhQUFGBoaAg0TSMnJ4crARYWFqK4uDihQCQ6iRZd+goEAlx1JNY5MqHNpiiqHeHQ+rssyx5Ldq6MhdnBYFCQAT2Bz+dDV1cXqqqqMD09zYXGQhAIBNDV1YXS0tJF9d+PP/6Ye7AkIhFRXjU1NYlaL+R3U61du1bU1c1qtXKadP41k5XKZDLBYrFApVJx4Xh0ArJ/fA4f9IzAw5sNzW/EEDOpo1EpsG1dFQJ2Q1KFW6ogiT+pVBqxNaJpGhaLBSaTievZJnttoVGfw+FAX18fGhsbI5Kg/L32l7/8Zfz85z9Pt7y46KagKOo8AL8GcIvQpouTotGC+IMRIcj09LTgY4n39rp16yKseUhYnZOTg46ODhQVFaG0tHTRGE4gHP4uLCwIVl4JBZGNZiJkX1hYwOTkJFpbWxdl2clKRYjodru5cJxhGBQXF6OkpAQ6nQ6bVldibWUJPugdQd/YHCfP5JsEioHNqyvRVFuMyfExtLa2iiopZVkWR48ehUKhQH19/aIuq1hJq+7uMD9IpSAnJyfmg9btdqO/v58beRRPsGKxWMQuTT0O4FKWZSeTvvI4Vnxlnp2dxdTUVMQf+OOPPxa0MhsMBk7HzP8io/fHNE3DZDLBaDTC7XZzjRC5ubkYGRkBTdNobGwUNbnj8Xg42Wj07Kqlguy9U51wCYT/LvzvIj8/nwvHZwwWPPXy26ApuWguk/k6DS5sb4BWFjZuaGlpEbXEx7Isl6xMNfIJBAIwm80wmUzc1oQk0WQyGdxuN3p6etDU1BSXqAzD4OOPP8bXvvY19PT0JPR6S4CIi6YoSg3ACWCC9+O3WZb9XwlPkiky0zSdsPuJ+Gf7fD5s3rw5ovSTjMykdm0ymdDS0pJSoothGFgsFiwsLGBhYQFqtRpr1qwR3AghBERkItT/SyhIuczj8WDz5s1LfvjwE0dmsxl+vx9lZeVwMAp0DE0vqTYtkVDYuqEWZzeuhtVixsTEBFpbW0WNfIjrqlarxZo1a5Z8Lrvdzm1NJBIJPB4PGhsbExK0o6MDd999N1566SXU1gp2+ImGKHuvFQmzg8Eguru7UVBQgA0bNsR8msbLyJI/oEQiWZQkE5LoIt5XLpcLGzduhEaj4ZRXpBc3XZkiEA5/JyYmRBeZkJZNmUzG1U2XChKOq9VqWCwWrF+/Pjxa1mhEU5kCIyYaRmcgZQKWFebiovaNKC3Iwfz8POfMKoZijIBhmIia+lLBH6ru8XjQ2dmJqqoqzMzMYHR0lBu2XlBQwN1zXV1duOuuu/D8888vhciiYdlX5nj2PHwcPHhw0dQIIBwWdXZ2xrTNFaroSiTYIDa0JpMJLMuiuLgYpaWlghI1pC3SbDajqalJ1Bs3FAqhp6cHBQUFgho8UoHL5UJvb+8iqWogEIDJZELP8Dg+6p9ACOFxLyqVKu77K2QyfG7zGrStq4ZEIsHMzAwWFhbQ0tIiquiGYRj09vYiLy8vbt03XZDt0aZNm7ioinRZmUwmWK1WdHZ2wmAw4OWXX8ZLL72E9esX+3GnCFH+oMtKZrLHTTZf6vDhw4vCZ6fTiZ6enpjzk4VMlQDC85snJyfR3NycNAFDbmaDwQCfz8d1OMUyHGBZFkNDQwiFQti4caPoTQ3d3d2orq4WbP4vFDabDYODg0kz+IFgEG8f7sOH3cfg8XhimunXVRThgq0NyNOGv9fJyUlYLBY0NzeLqhhjGIabipKo/TUdEIltY2NjXA0+y7J4/fXX8eMf/5gz73/66aeXGuaf3GQmOljghLTObDYvImksdHZ2oqGhgSOc0ERXImnm+Pg47HY7mpqaUl4liOGA0WiEw+GIqOECQG9vL3JycrBmzRpRV02ySohd9wZOlOJiZcPjwer04J9HBjEyrYfH44HH44FSLsUX2jbg7JYG7oEwPj4Ol8slyr6eDxKhkL5kMUGIvHHjxoSNQMPDw/jqV7+Kp59+Gk1NTbDb7TH7slPEqUFmvj67oaFB0B+3u7sba9euhVarjfsQSMVM4OjRo5DJZFi/fv2Sby4iqTQajTCZTAgEAigtLcW6detEDa0dDgf6+/tF7dQi0Ov1mJ6eFvRgjYXByfA4ndryQnyusRZuZ9iAwePxgKIoKBQKtLS0iK7h7u7uRmlpqehlPqJxSEbk8fFxfOlLX8KTTz6JtrY2MS/h5Cezy+VCV1cX58skFH19fZx4RCaTLXoICFV0BQIB9PT0oLS0VPSQjOw1a2pqEAwGYTQaIZVKuQTaUuqoZrOZG8ezVDP1aJCy1lL3saEQA6k08m8yMDAAv98PlUoFm80GrVbLiVWW8qALhULo6upCeXl53D74dEGI3NDQkLCuPjU1hRtvvBG///3v01I3JsHJTWaSKEjHCL+vrw92ux01NTWLSMgncqJVlsx5ykSdl4yeia4/+nw+GI1GGI1GBINBTpwRT5AQCyT7m+6qGQ9kq+F0OtHU1CRq+EuiH6VSyYk2iECDlL0oiuK+j1SUXzRNo7u7G5WVlaLnDIQSeXZ2Ftdffz0effTRlJSJKeDkJrPP54PP50t5hXI4HOjo6EBtbe0iaZzQRBdZ2TZv3ix6wzjpb25ubk641+QLVVyusIF9SUlJRGkjGhMTE7BarWnt6xOBmDuQBJ2Y+3qSWc7NzU1YIiJdTmQwAfk+8vPz434fNE1zkV2qLbDJ4Pf70dnZiQ0bNiRcbPR6Pa677jr8/Oc/x65du0S9Bh5ObjIzDCPIZZOPhYUFrqZXUFDA/QGF7o+BE6Z4zc3NoowkISBCFZvNljLZGIaB1WqF0WiE1WqFTqfjwk+ZTMaRLRgMiq5Ei7VqigWSkCoqKkppG0NKPUajETabDTqdjtNLk3CcjApatWqV6I0eQolsMBhw7bXX4uGHH8YFF1wg6jVE4eQmM8uyguc8kWkWVqsVLS0tmJ2dhVwuR1VVVUTPaaKpi6QzyefzYdOmTaKXQwYHBwFAcBIvHliWhdPp5BJocrkcwWAQeXl5cQU06YJfnxa7HkvC36XuY/nfh9ls5oQsCwsLWLt2bdp+a/FAtArr1q1L2ERiMplw7bXX4oc//CEuueQSUa8hBk4PMhMjAH6ia3p6GizLcraoyRJdoVAIvb290Ol0oncmkWkb+fn5i6yHxDh3Z2cnFAoFl/knCbRY7hmpgKjsKioqRE8akVWzpqZG9PDX6XSiu7sbcrkcLMtGDBRYasRCiFxfX5+w1Ge1WnHNNdfgwQcfxJVXXrmk9xSIU1fOSeD3+9HV1bXI/1oqlcLv9wsiss/ny8icJ/65V61alZE9G5nkQcJIkhUfHR2F1+vlbuRYg+qEnLuurk70lY2ce/Xq1ek2FcRFIBDA0aNH0dDQgOLiYoRCIVgsloj5Xsk8wBKdWwiR7XY7rr/+etx///3LRWTRkLGVGQj/4ePB4XCgt7cXDQ0Ni77c+fl5zM7Oor6+PmbLIv8c/f39aGhoEH1CIhkRk4lzk0z7+vXr44Z65EY2Go2w2+3Izc3lJmMk2kIQ8UOic8eC0wlMTFBgGKC2lkWs5C7J/mZCxEIe7PHIFh2OkzJgcXFx0vId6Xdfs2ZNwsqG0+nEddddh29961u48cYbl/yZUsDJHWYD8cms1+s5Hyx+mYLsj2mahtFohMFggNfr5XqR+e6UxEsrE7VYkg0X26gAOKENT6WjimVZ2O127kYmZgMlJSWLJK99fX0pC00cDuCll6QguyKJBNi9OwQ+p4gaLVkZJx0QIifbx/JByoAmkwl+vz+u3DYYDKKzszMpkd1uN2644QbcdtttuOWWW5b8mVLEyU/mQCAQMXArOtHFD5XiZayJlNJgMMDpdCI/P58zTk+nnzcZ5ubmMDs7K3o2HAgbDZIH0FJEJcRswGg0cvVblUqFiYmJtB5Ahw9T6O2VgOwkTKbw6rxrV7gfPV4zhhggq32yzHIiRMttyXyv3NxcbnBBoi2B1+vFjTfeiC996Uu4/fbb0/0oS8GpRWYi61QoFNiwYUNaii4yB9nn84GiKOTm5qK0tDSm7WyqIA8aIqoQMxsOhIUH8/Pzix5iS4Xf78f4+Djm5+cjVmwhHtsEn3xC4dgxCcjCZbMBpaUsLryQ4bYyiRr00wXZEoi52hNL3oWFBczMzECj0aCysjKuKs/n8+Hmm2/Gnj17cMcdd4hrWigcpw6ZyaC3ysrKRQJ5oUQmZvjFxcVcssxut8NgMMBisUCj0aC0tJSr3aYCvn5b7PIQUV45HI6MPCT4ijGJRMKtUCSKIQ0hiTLBej3wt79JkZMTDrFtNgqXXhpCbq4Vg4ODGdnKCG1sSAekSlBbW4ucnBxOrELmexFzv2AwiFtuuQUXXXQR7rrrrowT+bLLLkNlZSV+//vfR//q5M9mUxQFm80Wd9Ab37M40c1G9mv86Q9AeAIECbtdLhcMBgMmJyehUChQWlq6aE8ZC/yHhNgN5sRkDgBaWlpEv1nInKq2tjbuAUZM4/kuIseOHUuoky4vB664gkFXVzgBdvbZDHQ6MwYHh1PqqhIKMnMsE2E7IfKqVau4eyV6vtfs7Cy+9KUvwWq1oq2tDV//+tczTuQ33ngDXV1doldc+Mjoyjw1NYXR0VG0trZGPNlTUXRZreHVIZU/PDEZIHtKYugWfVN6vV709PSgrq5OdJVRrOkSYoFsCVwulyCdNV8nbTKZuExwrO8ECO/tx8fHRbf5AcD5aoltqQSckH/W1NQk/HvSNI1vfOMbKC4uRn5+Po4cOYLXX389Y4R2u9244IILcNNNN6Gvry9jK3NGyTwxMYGSkpKIsDcVIs/NzWFmZiapDjoRSNbTYDBwoozS0lJOrJKOUX8yEMFGoukS6YIYITAMk7bOmt8QQtN0hHPpwsICpqen0draKnpykSTSMknkZDruUCiEb37zm1i7di0eeuihZdkjf+Mb38AXv/hFeDwefPjhh6dmmF1dXR3hNpLqnCe3242tW7cuaZ+pUqm4MIuIMvr6+uByuVBZWQmJRBLXbywdkPzAmjVrRBdVEP8zlUq1pL199HdiMpkwPj4Om80GANi4caPoe3tC5Ewk0kivc1VVVVIi33333aiurl42Ij/55JOgKAo33HADnnjiiYy+V0ZXZr51kNBEF1kx1Wq16I0BQLgRQ6/XY9OmTXA4HDAYDHC5XCgsLERpaWnKais+yA2bidU+kzprILwlMplMqK6uhtlshs1m40o8S3UuJfVvsY3vgRO9zslaJBmGwbe//W3odDo88sgjojazJEJ7eztsNhtkMhnsdju8Xi+uueYa/OEPf+C/7OQPswmZhZrtkTlPlZWVouuJWTbs2+z1ehc1YhD7XYPBALvdjry8PK7kJfSPHm+6hBggYXtlZWVGEij8bDv5vKTEQ4Qq6TqXktJWS0uL6BlxsiKXl5cn/F4YhsEDDzwAAPif//mfZSNyNJ544olTN8wGhPcgO51O9Pf3pyxDFAJiy6pSqWLa1PLHgxJbIIPBgGPHjkGn03Elr3ihZ6LpEktFJnXW5AFHxp7yb3KKopCXl4e8vDzU19dHDKljWTaiISQe7HY7BgYGMkrksrKypET+3//7fyMQCODxxx9fMSIvBzK6Mj/55JNYs2YNWltbE+7BSHNBJle1srKylE3giB6Y2O+qVCqu5EWSQ0uZLpEMpCSXiQccSaSxLIuGhoaUJ0Ekcy4lzp/pjK5NBj6RE0VwLMviRz/6Eebn57F3717R8wAi4uQPs1944QX85S9/wdDQEM4//3xcddVV2LZtW0QoR2qlTU1NopdBCBnWrl0rSjKKlHeI3xeJNohgQ0ykq7MWgkSzmVJFLOdStVoNvV4v+iAAILzSdnd3o6SkJGGlgGVZPPzwwxgZGcGTTz4pqnNLBnDyk5nA6/Xi9ddfx/79+9Hd3Y1du3bh8ssvxyuvvIKbbroJW7ZsEZ0MiczulwpilRMIBDi/K1LyEiOcFOpnnQ7IlkOn02Wk/j09PY2xsTEoFApOqMKPZJYCQuRkVrssy+KXv/wlOjs78Ze//EX0iCkDOHXIzIff78eLL76Ie++9F6WlpWhra8M111yDnTt3ivalGwwGjI+PL7mhIRZiZZUDgQAMBgMMBgNn5FdWVpaWwYDJZMLo6KjoA9b4115YWJiRcSpms5nz4lYoFBENIUt1LiXm90VFRUmJ/Nhjj+HDDz/Ec889J3q0lyGcmmQGgB/84AdoaWnBlVdeiXfeeQcHDhzARx99hO3bt2PPnj3YtWtX2n+EqakpGI3GjOxhhUyXIHXbRO2b8ZApZ07ghM1PWVmZ6EIW4MRDqK2tLea1+3w+7ntJ1bmUREIFBQUJvcZYlsXevXvx5ptv4sCBA6J3vWUQpy6ZY4GmaXzwwQfYt28f3nvvPbS1tWHPnj04//zzBa1QxBQvEAhg06ZNooft6UyXiG7fLCgo4GrZ0ddHcgfNzc2i7++IzU8mRtwAqcs/iUaa1PgTOZcSIufn5yeNJp544gm89NJLeOmll0SPajKM04vMfIRCIXz88cfYv38/3n77bTQ2NmLPnj246KKLYu5JiQ5aq9WK7gEGiDNdgjh0GgwG2Gw2rn2zoKAAk5OTcLvdoo9zAU64bGSitAWEtzRk6mU6kVAi51KJRIK+vj7k5uYmFco8/fTTeOaZZ/Dyyy+LXgZbBpy+ZOaDYRgcPnwY+/btwz/+8Q/U19dj9+7duOSSS5CTkwOn04nBwcGMCE2AzEyXIM4hBoMBc3NzkMlkXMZdzFWZSEuT+V6li4WFBUxNTYmm4452LvX5fMjLy0NDQ0PClXbfvn34wx/+gFdffVV0qegy4cwgMx8Mw6Crqwv79+/H3//+dxQVFWF8fBx/+tOfxJ79AyCze1i+zrqsrIy7gVNp30wE0i+8FAePRNDr9ZiZmUFra6vo2wKWZdHX1welUgmVSgWj0RjXufTFF1/EY489hldeeUX0qgVBIBDAPffcg7feegssy+KnP/0prr32WjHf4swjMx8fffQRbr/9dpx//vk4dOgQiouLsWfPHlx++eWirEKZmi4BJM4qC23fTATSZpiJGjVwwnAxU0Tu7++HRqOJGJManVh87733IJPJ8Nprr+G1117LyAOLQK/X48MPP8R1112H4eFhbN++HUajUcwEa5bMtbW1qK6u5tRM+/fvx8svv4zc3Fzs3r0bV155JUpKSlLaQ2dyugSQms46XvtmotozEZtkojsJCLelzs/PJ1X1pQMiZlGpVItGE/ERCoXws5/9DE8//TQUCgW2b9+O3/3ud8tWTy4uLsbo6KiYkcCZTeZ4IO2TBw4cwEsvvQSFQoHdu3fjqquuQnl5eUJiE0GFRqPJSCKNuFCuXr065WQUad8kEsri4mKUlpZGlHaI2CQTNj9AuOOMTJDMFJGVSmXS7/6dd97BQw89hFdffRXFxcXo6urCli1bRL2eePjjH/+Ip556Cm+//baYp82SORlYlsXU1BQOHDiAF198EQzD4Morr8SePXtQXV0dccOQOmxpaanog7wBcXXW0aWdwsJCqNVqLvQVWygDhEtnJpMJzc3NGSHywMAA5HJ5UnnpBx98gO9973t49dVXRR9MkAw//elP8eyzz+K1114Tu8SXJXMqYFkW8/PzOHDgAF544QV4vV5cfvnluOqqqyCVSnH48GGcd955otsHAZnVWTMMg/HxcUxPT0OhUCA/Pz/l9s1kmJqagsViQXNzs+jbDuKTJpPJkhL5k08+wb333otXXnklI5WLRPi3f/s3uN1uPProo5mIerJkXgoMBgNeeOEFPPXUUxgeHsZ1112Hf/mXf8H69etFDa9tNhsGBgYy0pgPnCgPEQtf0r5psVgEtW8mw+TkJDf5MlNElkqlWLduXcLvvaOjA3fffTf+9re/pTRxUgx8+umnePDBB/HWW29l6i2yZF4qZmdncfnll+NXv/oVjh07hgMHDkCv1+Piiy/G1VdfjY0bNy7pBiam95nQWQPhZNTc3FzMrLKQ9s1kIMPZMyFmIUlLiqKSPkC7urrwzW9+Ey+88EJEhnu58Pjjj+OBBx6ImIjx61//WszpkFkyLxUsy8JqtUbsYW02G15++WU8//zzGB8fx0UXXYQ9e/ak3OaYyRo1cEL+KTQZxW/flMlkXGY8nn55bGwMbrc7I9JYUjFgWTapl1lfXx++8Y1vYP/+/Vi/fr2o13ESIUvmTMPpdOLVV1/FgQMHMDQ0hAsuuABXXXUV2tvbE97gxE8rEzprIFwDt9vtaYe+Xq+Xy4xHt28SG19iryR2Rp/M0WYYJimRBwYG8PWvfx3PPPMMGhsbRb2OkwxZMi8n+D3ZPT092LVrF6666iqcffbZ3MpIiJApnTUpuxGiiXH+6PZNiUQCuVyesWTXyMgIaJpO6m4yPDyMr371q3j66afR1NQk6nWchMiSeaXg8/nwj3/8A/v378eRI0ewY8cO7N69G6+++ipuvPFGbNu2LSMr2vDwMEKhUNp+2cnOPzQ0BJfLBblcnnL7ppDzj4yMIBgMJr3+8fFx3HzzzXjiiScyItM9CXFqkPlXv/oVHn/8cW6I9X//93+v1HCujCAQCOAf//gH7rnnHmg0GmzZsgVXX301zjvvPNH2yqQOK5VKRc+2k/NH72FTad8Ucv7R0VH4/X40NjYmvP6pqSnceOON+P3vf49t27Yt9aOdKjg1yDw3N4fKykp4vV40Njbi5ZdfxubNm5d62pMKv/zlL8GyLL71rW/h/fffx759+/DBBx9wPdlf+MIX0s5mk4YMolXOBJEHBwchkUjiPijitW8Knb45OjoKn8+XlMizs7O44YYb8Jvf/AY7duxY0uc6xXBqkJlgcnISl156KT7++GPRDeJXGrEmYoRCIXz00Uc4cOAA3n77bWzatAl79uzBhRdeKFh0EAqFuMb8TBjfkxVfJpMlrfPyjyHtm2azGVqtNuH0zbGxMXg8nqTJNL1ej+uuuw6/+MUvcN555y3pc52CODXIfOzYMVx00UUwGAx49NFHceutty71lKccGIbBoUOHsH//fq4ne8+ePbj44ovjNkMQO9mSkpKMyEtT0UInOgeZvhmrfXN8fBwulwubN29OeH6DwYBrr70WDz/8MC644IKlfKxTFScXme+44w4cOXIk4md79+5FS0sLgHBd9PLLL8ejjz6Kc845J51rPS1AerL37duH119/HatWrcLu3btx2WWXcV04pLOqqqoqIzY/pM1QrVaLGrrz2zf9fj9kMlnSpg+TyYRrr70WP/zhD8UUYcTEc889h/vvvx9SqRTf+973cNttt2X0/VLAyUVmIXjwwQeRn5+Pe++9V8zTnrIgTfj79u3Da6+9hpKSElx44YV47bXX8Nvf/jYjo2j4VruZUlORXvCioqIIY4Ho9k2r1YprrrkGDz74IK688sqMXAuB0+lEY2MjPv30U0ilUrS2tqK3t1f04X5p4tQg80cffYSdO3fC5XLhggsuwH/913/h85//fErneOGFF/CDH/wATqcTF1xwAX73u9+dzNMJ0gLLsvjggw9w8803Y9WqVdDpdNi9ezeuuOKKlHuy44EQmcyMzgRiabmjJ2AYjUaoVCr89Kc/xXe/+11cc801GbkWPvbv348XX3wRf/7znwEAN998M3bv3o2bbrop4+8tAKKQOeODd37yk59g1apV2Lp1K2688caUiQyEn+AHDx7EsWPHMDk5iWeffVb8C11hUBSFTz/9FH/5y1/w0Ucf4dFHH4Xb7cbNN9+MK664Ao8//jjm5+eR5OEbF8TlMi8vL2NEnpqaitmUoVAoUFlZidbWVrS3tyMQCOD73/8+5ubm8OGHH2JsbCwj18PH9PR0hKtLdXU15ufnM/6+y4mMz+x49dVXl3wO/t6mpaUFRqNxyec8GXHfffdx/11fX4/7778f9913H9eTTZKHV1xxRcye7HggkyCKiooy1nEktE3S7/fjD3/4A+6//37ccMMNYjf5x0UgEIi4LolEctpFd6fUSLy5uTk8//zzuOyyy1b6UpYNFEWhtrYW3/72t/H+++/j2WefhUajwZ133okLL7wQP//5zzE2NhZ3xSZZ8eLi4owReXp6GmazOSmRvV4vbrrpJtxyyy245ZZboFQqcemlly5LJ1RFRQVmZ2e5/5+ZmclIlWAlcVLJORNlxDs6OnDzzTfjJz/5Ca677rrlvKyTEizLcj3Zzz//PGw2Gy677DJcddVVnPiDELm0tDQjUyyAMClI91YiIvt8Ptx8883Ys2cP7rjjjmVXAS4sLGDLli3o7OwEwzDYsWMHent7M9JjngZOjQSYGHj//ffxrW99C0899RRX6soiEmazGS+++CKef/55LCws4Pzzz8cHH3yAX/ziFxlrVBDqCRYIBPCVr3wFX/ziF3HXXXetmJz3iSeewA9/+EMAwCOPPIKrr756Ra4jBs4cMjc2NuLvf/+7KMPODh48iH/7t3/Diy++mLHVaqUxPT2Niy++GCUlJXA4HFxPtpidULOzs1hYWEhK5GAwiFtvvRU7d+7Ed77zndNKly8izgwye71e5OTkRMgZzz//fPzud79L+Vz33Xcfurq60NnZic7OztOWzC+//DI8Hg9uvPHGiJ7s4eFhrid769ataRN7bm4Oer0+KZFpmsbtt9+OtrY2/Md//EeWyPFxZpBZTNhsNk7n/OGHH562ZI4Hj8eDv//97zhw4AD6+vq4nuyzzjpLcGZ3fn6esypKdEwoFMKdd96J+vp6PPTQQ1kiJ0aWzOniTCUzH/ye7M8++ww7duzA1VdfjR07dsR1R0mFyHfffTfKy8vxk5/8JEvk5MiSOR6S6cSzZI5EIBDA22+/jQMHDuCTTz7BWWedhT179uDcc8/lerKFzpZiGAbf/va3odPp8Mgjj4juVnKaIkvmdJElc3zQNB3Rk71lyxaUlZXB6XTi4YcfTkrkBx54AADwP//zP1kiC8epIec8XfDcc89h9erVqK+vxx/+8IeVvpyMQSaT4fzzz8djjz2G7u5urF+/Hs888wwOHjyIO+64A3/729/g8XgWHccwDH7wgx8gEAhkibxCyLic83SA0+nEd77znYiOGzKU7nQGMSgkzRmHDh3Cvn378J//+Z9Yt24d9uzZgy9+8YvQarX40Y9+BIvFgr1792aJvEI4I8PsVHGSd9wsOxiGQWdnJ/bt24c33ngDgUAA69evx/79+087vfMyQZQwO7syC8CZ0HGTCiQSCbZu3YqtW7fiJz/5CV555RWcf/75WSKvMLLxkACcCR036UIikWD37t0ZmQUdCy+88AKamppQV1eH22+/HaFQaFne91RAlswCcCZ03JwqOBN629NFlswCcPHFF+ONN96AwWCAXq/Hxx9/jC9+8Ytpncvv9+Oxxx47mUT+pxRuu+02aDQayOXy07q3PR1k98wCUFZWhh//+Mf43Oc+BwD42c9+lnbr3IYNG9DW1gan0ynmJZ5xIL3tb7755kpfykmDbDZ7mWGz2dDV1YUf/ehHmZz3e8rjDOttz2azT0WcbgMAMoXf/va3MX9Oetv37duX7W2PQpbMWZxSuPPOO0XrbT/dkCVzFqcMvF4v15NNkG5v++mILJlPQQQCAdxzzz146623wLIsfvrTn+Laa69d6cvKONRqNWiaXunLOGmRLU2dgrBYLDj//PMxPDyMV199FbfffjuCweBKX1YWK4xsNvs0QHFxMUZHR7lZVVmccsi2QGYB/PGPf0Rzc3OWyFlk98ynMn7605/i2WefxWuvvbbSl5LFSYBkYXYWJykoivoNAC2Af2VZdrFbQBZnHLJh9ikIiqLOBrCBZdlbl0pkiqIkFEX9g6KoYYqihiiKuliky8ximZEl86mJVgDtFEWN8P5Jd1I5C+CrLMuuB/DvAH4s1kVmsbzIhtlZcKAo6g4A21mWvX2lryWL1JFNgGUBiqLuA3A/ACOAbJh9iiK7MmfBgaKoawD8BMBGNntjnHLIkjmLCFAUNQOglWVZ00pfSxapIZsAO8NBUdQaiqLKj//35wD4skQ+NZHdM2eRD+B1iqKkAAwAblzZy8kiXWTD7CyyOE2QDbOzyOI0QZbMWWRxmiBL5iyyOE2QJXMWWZwmyJI5iyxOE2TJnEUWpwmyZM4ii9MEWTJnkcVpgiyZs8jiNMH/D0YUjDRB3n3nAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "迭代次数：2095次\n",
      "梯度：[ 7.33881397e-05  2.44073067e-05  2.52604176e-04 -5.13424350e-04]\n",
      "权重：[  4.34496173   2.05340452   9.64074166 -22.85079478]\n",
      "模型准确率：100.00%\n"
     ]
    }
   ],
   "source": [
    "# 训练数据集\n",
    "X_train = np.array([[3, 3, 3], [4, 3, 2], [2, 1, 2], [1, 1, 1], [-1, 0, 1], [2, -2, 1]])\n",
    "y_train = np.array([[1, 1, 1, 0, 0, 0]])\n",
    "# 构建实例，进行训练\n",
    "clf = MyLogisticRegression(epsilon=1e-6)\n",
    "clf.fit(X_train, y_train)\n",
    "clf.draw(X_train, y_train)\n",
    "print(\"迭代次数：{}次\".format(clf.n_iter_))\n",
    "print(\"梯度：{}\".format(clf.grad_))\n",
    "print(\"权重：{}\".format(clf.w[0]))\n",
    "print(\"模型准确率：%0.2f%%\" % (clf.score(X_train, y_train) * 100))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 习题6.3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "&emsp;&emsp;写出最大熵模型学习的DFP算法。（关于一般的DFP算法参见附录B）"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**解答：**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**解答思路：**\n",
    "1. 写出最大熵模型\n",
    "2. 根据附录B的DFP算法，写出最大熵模型学习的DFP算法\n",
    "3. 自编程实现最大熵模型学习的DFP算法"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**解答步骤：**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**第1步：最大熵模型**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "&emsp;&emsp;根据书中第6.6.2节的最大熵模型的定义：\n",
    "\n",
    "> **定义6.3（最大熵模型）** 假设满足所有约束条件的模型集合为\n",
    "> $$ \\mathcal{C} = \\{P \\in \\mathcal{P} | E_P(f_i) = E_{\\tilde{P}}(f_i), \\quad i=1,2,\\cdots,n \\} $$\n",
    "> 定义在条件概率分布$P(Y|X)$上的条件熵为\n",
    "> $$ H(P) = -\\sum \\limits_{x,y} \\tilde{P}(x) P(y|x) \\log P(y|x) $$\n",
    "> 则模型集合$\\mathcal{C}$中条件熵$H(P)$最大的模型称为最大熵模型。式中的对数为自然对数。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "&emsp;&emsp;根据书中第6.2.2节的最大熵模型的特征函数：\n",
    "\n",
    "> 用特征函数$f(x,y)$描述输入$x$与输出$y$之间的某一个事实。其定义是\n",
    "> $$ f(x, y) = \\left \\{ \\begin{array}{ll}\n",
    "1, \\quad x与y满足某一事实 \\\\\n",
    "0, \\quad 否则\n",
    "\\end{array} \\right. $$\n",
    "它是一个二值函数，当$x$和$y$满足这个事实时取值为1，否则取值为0。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**第2步：最大熵模型学习的DFP算法**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "&emsp;&emsp;根据书中第6.3.2节拟牛顿法的最大熵模型学习：\n",
    "> 对于最大熵模型而言，\n",
    "> $$ P_w(y | x)=\\frac{\\displaystyle \\exp \\left(\\sum_{i=1}^n w_i f_i (x, y)\\right)}{\\displaystyle \\sum_y \\exp \\left(\\sum_{i=1}^n w_i f_i(x, y)\\right)} $$\n",
    "> 目标函数：\n",
    "> $$ \\min \\limits_{w \\in R^n} \\quad f(w) = \\sum_{x} \\tilde{P}(x) \\log \\sum_{y} \\exp \\left(\\sum_{i=1}^{n} w_{i} f_{i}(x, y)\\right)-\\sum_{x, y} \\tilde{P}(x, y) \\sum_{i=1}^{n} w_{i} f_{i}(x, y) $$\n",
    "> 梯度：\n",
    "> $$ g(w) = \\left( \\frac{\\partial f(w)}{\\partial w_1}, \\frac{\\partial f(w)}{\\partial w_2}, \\cdots, \\frac{\\partial f(w)}{\\partial w_n}\\right)^T $$\n",
    "> 其中\n",
    "> $$ \\frac{\\partial f(w)}{\\partial w_i} = \\sum \\limits_{x,y} \\tilde{P}(x) P_w(y|x) f_i(x,y) - E_{\\tilde{P}}(f_i), \\quad i=1,2,\\cdots,n $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "&emsp;&emsp;根据书中附录B的DFP算法：\n",
    "\n",
    "> **算法B.2（DFP算法）**  \n",
    "> 输入：目标函数$f(x)$，梯度$g(x) = \\nabla f(x)$，精度要求$\\varepsilon$；  \n",
    "输出：$f(x)$的极小值点$x^*$  \n",
    "（1）选定初始点$x^{(0)}$，取$G_0$为正定对称矩阵，置$k=0$。  \n",
    "（2）计算$g_k=g(x^{(k)})$，若$\\|g_k\\| < \\varepsilon$，则停止计算，得近似解$x^*=x^{(k)}$，否则转步骤（3）。    \n",
    "（3）置$p_k=-G_k g_k$。  \n",
    "（4）一维搜索：求$\\lambda_k$使得\n",
    "> $$\n",
    "f (x^{(k)}+\\lambda_k p_k )=\\min \\limits_{\\lambda \\geqslant 0} f(x^{(k)}+\\lambda p_{k})\n",
    "$$\n",
    ">（5）置$x^{(k+1)}=x^{(k)}+\\lambda_k p_k$。  \n",
    "（6）计算$g_{k+1}=g(x^{(k+1)})$，若$\\|g_{k+1}\\| < \\varepsilon$，则停止计算，得近似解$x^*=x^{(k+1)}$；否则，按照式（B.24）算出$G_{k+1}$。  \n",
    "（7）置$k=k+1$，转步骤（3）。  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "&emsp;&emsp;根据公式（B.24），DFP算法的$G_{k+1}$的迭代公式为：\n",
    "\n",
    "$$\n",
    "G_{k+1}=G_k+\\frac{\\delta_k \\delta_k^T}{\\delta_k^T y_k}-\\frac{G_k y_k y_k^T G_k}{y_k^T G_k y_k}\n",
    "$$\n",
    "\n",
    "&emsp;&emsp;其中$y_k= g_{k+1} - g_k, \\delta_k = w^{(k+1)} - w^{(k)}$。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "最大熵模型的DFP算法：   \n",
    "输入：目标函数$f(w)$，梯度$g(w) = \\nabla f(w)$，精度要求$\\varepsilon$；  \n",
    "输出：$f(w)$的极小值点$w^*$  \n",
    "（1）选定初始点$w^{(0)}$，取$G_0$为正定对称矩阵，置$k=0$。  \n",
    "（2）计算$g_k=g(w^{(k)})$，若$\\|g_k\\| < \\varepsilon$，则停止计算，得近似解$w^*=w^{(k)}$，否则转步骤（3）。  \n",
    "（3）置$p_k=-G_k g_k$。  \n",
    "（4）一维搜索：求$\\lambda_k$使得\n",
    "$$\n",
    "f(w^{(k)}+\\lambda_k p_k) = \\min \\limits_{\\lambda \\geqslant 0} f(w^{(k)}+\\lambda p_k)\n",
    "$$\n",
    "（5）置$w^{(k+1)}=w^{(k)}+\\lambda_k p_k$。  \n",
    "（6）计算$g_{k+1}=g(w^{(k+1)})$，若$\\|g_{k+1}\\| < \\varepsilon$，则停止计算，得近似解$w^*=w^{(k+1)}$；否则，按照DFP算法的迭代公式算出$G_{k+1}$。  \n",
    "（7）置$k=k+1$，转步骤（3）。  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**第3步：自编程实现最大熵模型学习的DFP算法**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "import copy\n",
    "from collections import defaultdict\n",
    "\n",
    "import numpy as np\n",
    "from scipy.optimize import fminbound\n",
    "\n",
    "\n",
    "class MaxEntDFP:\n",
    "    def __init__(self, epsilon, max_iter=1000, distance=0.01):\n",
    "        \"\"\"\n",
    "        最大熵的DFP算法\n",
    "        :param epsilon: 迭代停止阈值\n",
    "        :param max_iter: 最大迭代次数\n",
    "        :param distance: 一维搜索的长度范围\n",
    "        \"\"\"\n",
    "        self.distance = distance\n",
    "        self.epsilon = epsilon\n",
    "        self.max_iter = max_iter\n",
    "        self.w = None\n",
    "        self._dataset_X = None\n",
    "        self._dataset_y = None\n",
    "        # 标签集合，相当去去重后的y\n",
    "        self._y = set()\n",
    "        # key为(x,y), value为对应的索引号ID\n",
    "        self._xyID = {}\n",
    "        # key为对应的索引号ID, value为(x,y)\n",
    "        self._IDxy = {}\n",
    "        # 经验分布p(x,y)\n",
    "        self._pxy_dic = defaultdict(int)\n",
    "        # 样本数\n",
    "        self._N = 0\n",
    "        # 特征键值(x,y)的个数\n",
    "        self._n = 0\n",
    "        # 实际迭代次数\n",
    "        self.n_iter_ = 0\n",
    "\n",
    "    # 初始化参数\n",
    "    def init_params(self, X, y):\n",
    "        self._dataset_X = copy.deepcopy(X)\n",
    "        self._dataset_y = copy.deepcopy(y)\n",
    "        self._N = X.shape[0]\n",
    "\n",
    "        for i in range(self._N):\n",
    "            xi, yi = X[i], y[i]\n",
    "            self._y.add(yi)\n",
    "            for _x in xi:\n",
    "                self._pxy_dic[(_x, yi)] += 1\n",
    "\n",
    "        self._n = len(self._pxy_dic)\n",
    "        # 初始化权重w\n",
    "        self.w = np.zeros(self._n)\n",
    "\n",
    "        for i, xy in enumerate(self._pxy_dic):\n",
    "            self._pxy_dic[xy] /= self._N\n",
    "            self._xyID[xy] = i\n",
    "            self._IDxy[i] = xy\n",
    "\n",
    "    def calc_zw(self, X, w):\n",
    "        \"\"\"书中第100页公式6.23，计算Zw(x)\"\"\"\n",
    "        zw = 0.0\n",
    "        for y in self._y:\n",
    "            zw += self.calc_ewf(X, y, w)\n",
    "        return zw\n",
    "\n",
    "    def calc_ewf(self, X, y, w):\n",
    "        \"\"\"书中第100页公式6.22，计算分子\"\"\"\n",
    "        sum_wf = self.calc_wf(X, y, w)\n",
    "        return np.exp(sum_wf)\n",
    "\n",
    "    def calc_wf(self, X, y, w):\n",
    "        sum_wf = 0.0\n",
    "        for x in X:\n",
    "            if (x, y) in self._pxy_dic:\n",
    "                sum_wf += w[self._xyID[(x, y)]]\n",
    "        return sum_wf\n",
    "\n",
    "    def calc_pw_yx(self, X, y, w):\n",
    "        \"\"\"计算Pw(y|x)\"\"\"\n",
    "        return self.calc_ewf(X, y, w) / self.calc_zw(X, w)\n",
    "\n",
    "    def calc_f(self, w):\n",
    "        \"\"\"计算f(w)\"\"\"\n",
    "        fw = 0.0\n",
    "        for i in range(self._n):\n",
    "            x, y = self._IDxy[i]\n",
    "            for dataset_X in self._dataset_X:\n",
    "                if x not in dataset_X:\n",
    "                    continue\n",
    "                fw += np.log(self.calc_zw(x, w)) - \\\n",
    "                    self._pxy_dic[(x, y)] * self.calc_wf(dataset_X, y, w)\n",
    "\n",
    "        return fw\n",
    "\n",
    "    # DFP算法\n",
    "    def fit(self, X, y):\n",
    "        self.init_params(X, y)\n",
    "\n",
    "        def calc_dfw(i, w):\n",
    "            \"\"\"计算书中第107页的拟牛顿法f(w)的偏导\"\"\"\n",
    "\n",
    "            def calc_ewp(i, w):\n",
    "                \"\"\"计算偏导左边的公式\"\"\"\n",
    "                ep = 0.0\n",
    "                x, y = self._IDxy[i]\n",
    "                for dataset_X in self._dataset_X:\n",
    "                    if x not in dataset_X:\n",
    "                        continue\n",
    "                    ep += self.calc_pw_yx(dataset_X, y, w) / self._N\n",
    "                return ep\n",
    "\n",
    "            def calc_ep(i):\n",
    "                \"\"\"计算关于经验分布P(x,y)的期望值\"\"\"\n",
    "                (x, y) = self._IDxy[i]\n",
    "                return self._pxy_dic[(x, y)]\n",
    "\n",
    "            return calc_ewp(i, w) - calc_ep(i)\n",
    "\n",
    "        # 算出g(w)，是n*1维矩阵\n",
    "        def calc_gw(w):\n",
    "            return np.array([[calc_dfw(i, w) for i in range(self._n)]]).T\n",
    "\n",
    "        # （1）初始正定对称矩阵，单位矩阵\n",
    "        Gk = np.array(np.eye(len(self.w), dtype=float))\n",
    "\n",
    "        # （2）计算g(w0)\n",
    "        w = self.w\n",
    "        gk = calc_gw(w)\n",
    "        # 判断gk的范数是否小于阈值\n",
    "        if np.linalg.norm(gk, ord=2) < self.epsilon:\n",
    "            self.w = w\n",
    "            return\n",
    "\n",
    "        k = 0\n",
    "        for _ in range(self.max_iter):\n",
    "            # （3）计算pk\n",
    "            pk = -Gk.dot(gk)\n",
    "\n",
    "            # 梯度方向的一维函数\n",
    "            def _f(x):\n",
    "                z = w + np.dot(x, pk).T[0]\n",
    "                return self.calc_f(z)\n",
    "\n",
    "            # （4）进行一维搜索，找到使得函数最小的lambda\n",
    "            _lambda = fminbound(_f, -self.distance, self.distance)\n",
    "\n",
    "            delta_k = _lambda * pk\n",
    "            # （5）更新权重\n",
    "            w += delta_k.T[0]\n",
    "\n",
    "            # （6）计算gk+1\n",
    "            gk1 = calc_gw(w)\n",
    "            # 判断gk1的范数是否小于阈值\n",
    "            if np.linalg.norm(gk1, ord=2) < self.epsilon:\n",
    "                self.w = w\n",
    "                break\n",
    "            # 根据DFP算法的迭代公式（附录B.24公式）计算Gk\n",
    "            yk = gk1 - gk\n",
    "            Pk = delta_k.dot(delta_k.T) / (delta_k.T.dot(yk))\n",
    "            Qk = Gk.dot(yk).dot(yk.T).dot(Gk) / (yk.T.dot(Gk).dot(yk)) * (-1)\n",
    "            Gk = Gk + Pk + Qk\n",
    "            gk = gk1\n",
    "\n",
    "            # （7）置k=k+1\n",
    "            k += 1\n",
    "\n",
    "        self.w = w\n",
    "        self.n_iter_ = k\n",
    "\n",
    "    def predict(self, x):\n",
    "        result = {}\n",
    "        for y in self._y:\n",
    "            prob = self.calc_pw_yx(x, y, self.w)\n",
    "            result[y] = prob\n",
    "\n",
    "        return result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "模型训练迭代次数：878次\n",
      "模型权重：[ 11.72083082  -0.7189626    9.69699232  -8.83693899   5.57382082\n",
      "  19.50768089   0.7189626   -9.69699232   8.83693899  -4.5237319\n",
      "   9.15646808  -6.6123125   12.96011049   4.5237319    6.6123125\n",
      " -12.96011049  -5.57382082  -9.15646808 -11.72083082]\n",
      "预测结果： {'no': 2.097720173362797e-16, 'yes': 0.9999999999999998}\n"
     ]
    }
   ],
   "source": [
    "# 训练数据集\n",
    "dataset = np.array([['no', 'sunny', 'hot', 'high', 'FALSE'],\n",
    "                    ['no', 'sunny', 'hot', 'high', 'TRUE'],\n",
    "                    ['yes', 'overcast', 'hot', 'high', 'FALSE'],\n",
    "                    ['yes', 'rainy', 'mild', 'high', 'FALSE'],\n",
    "                    ['yes', 'rainy', 'cool', 'normal', 'FALSE'],\n",
    "                    ['no', 'rainy', 'cool', 'normal', 'TRUE'],\n",
    "                    ['yes', 'overcast', 'cool', 'normal', 'TRUE'],\n",
    "                    ['no', 'sunny', 'mild', 'high', 'FALSE'],\n",
    "                    ['yes', 'sunny', 'cool', 'normal', 'FALSE'],\n",
    "                    ['yes', 'rainy', 'mild', 'normal', 'FALSE'],\n",
    "                    ['yes', 'sunny', 'mild', 'normal', 'TRUE'],\n",
    "                    ['yes', 'overcast', 'mild', 'high', 'TRUE'],\n",
    "                    ['yes', 'overcast', 'hot', 'normal', 'FALSE'],\n",
    "                    ['no', 'rainy', 'mild', 'high', 'TRUE']])\n",
    "\n",
    "X_train = dataset[:, 1:]\n",
    "y_train = dataset[:, 0]\n",
    "\n",
    "mae = MaxEntDFP(epsilon=1e-4, max_iter=1000, distance=0.01)\n",
    "# 训练模型\n",
    "mae.fit(X_train, y_train)\n",
    "print(\"模型训练迭代次数：{}次\".format(mae.n_iter_))\n",
    "print(\"模型权重：{}\".format(mae.w))\n",
    "\n",
    "result = mae.predict(['overcast', 'mild', 'high', 'FALSE'])\n",
    "print(\"预测结果：\", result)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 参考文献"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "【1】指数分布族（来源于Wiki百科）：https://en.wikipedia.org/wiki/Exponential_family  \n",
    "【2】逻辑斯谛不属于指数分布族的证明方法：https://stats.stackexchange.com/questions/275773/does-logistic-distribution-belongs-to-exponential-family  \n",
    "【3】逻辑斯谛分布（来源于Wiki百科）：https://en.wikipedia.org/wiki/Logistic_distribution  \n",
    "【4】指数倾斜（来源于Wiki百科）：https://en.wikipedia.org/wiki/Exponential_tilting  "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.5"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
